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FOUR METHODS OF ATTITUDE DETERMINATION 
FOR SPIN-STABILIZED SPACECRAFT 
WITH APPLICATIONS AND COMPARATIVE RESULTS 


Gene A. Smith 

Goddard Space Flight Center 

FORMULATION OF PROBLEM 
Spacecraft Attitude 

The attitude of a spacecraft is defined by the specification of three rotation angles relating 
the spacecraft axes to a known coordinate system. These rotation angles are Eulerian, and 
as such are not unique but are chosen for convenience or naturalness. (For further discus- 
sion of Eulerian angles and their nonuniqueness, see Appendix A and references 1 and 2.) 

A number of well defined coordinate systems are available for attitude design considera- 
tions. The most convenient of these systems, are those with inertial space coordinates, which 
have the advantage of being constant with time and are related to the celestial sphere. The 
two most common systems are the equatorial inertial (Z = north pole, X = first point of 
Aries, and Y = perpendicular to XZ plane to form a right-handed system) and the ecliptic 
inertial (Z = ecliptic north pole, X = first point of Aries, and Y as above). 

The spacecraft axes are chosen with regard to the design, function, moments of inertia, and 
planned orientation of the spacecraft, and are fixed with respect to the satellite body. As 
examples, a pointing spacecraft has the pointing direction defined as one axis (usually the 
X-axis), and the spinning spacecraft, discussed in this paper, has its spin axis defined as the 
Z-axis. 

Thus, with selection of the body-fixed spacecraft axes, the appropriate external coordinate 
system, and a set of rotation angles, the spacecraft attitude can be described unambiguously. 

Spin-stabilization 

If a spacecraft is spun about its largest principal moment of inertia, this spin axis will point 
in the same direction in inertial space unless acted upon by torque-producing forces. For 
most earth-orbiting spacecraft, gravitational forces are nearly constant, with the slight 
deviations adding to other variant forces to create spacecraft torques. These are relatively 
small perturbations which give rise to slowly varying spin-axis orientations. By slowly 
varying is meant a time much longer than the spin period, such that the spacecraft orien- 
tation does not become unstable. This condition is called spin stabilization. 

The design of the satellite and the proposed orbit are critical considerations in preventing 
torques from becoming disruptive to the spin stability. Thus, antennas and booms are kept 



light and short, to keep gravitational torques small; magnetic moments are balanced; and 
the orbit is high enough to make air drag negligible. In addition, some error in alignment 
of the maximum moment of inertia with the symmetric body axis, about which the space- 
craft is spun, is common and produces a coning of the spin axis about the net angular 
momentum vector. (See reference 2, pp. 113-125.) In this report it is assumed that exter- 
nal torques are negligible, that initial spacecraft oscillations have been damped, and that 
misalignment effects are less than other error sources. 

Spin-stabilized Spacecraft Attitude 

In the case of a spacecraft spinning about its largest principal moment of inertia, the spin 
axis is defined as the z-axis (figure 1). To transform the orthogonal spacecraft coordinates, 
xyz, into the orthogonal equatorial inertial coordinates, XYZ, choose the Eulerian rota- 
tion angles 0, 5 , and -a where they are defined as follows: First, rotate the x-axis counter- 
clockwise (ccw) by the arc 0 around the z-axis, to bring the x-axis into the Zz plane and 
the y-axis into the XY plane and perpendicular to the Zz plane. Second, rotate ccw around 
the new y-axis by b to make the z-axis coincident with the Z. Third, rotate clockwise (cw) 
around the final z (=Z) axis by -a to complete the transformation. 


Z NORTH 



Figure 1 . Spacecraft Axes, x, y, z, Referenced to Inertial Coordinates, X, Y, Z 
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In matrix notation, the transformation is 
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R(t) = 
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For the transformation of the z-axis only (that is, X, Y, and Z components of a vector 
parallel to the z-axis and x = y = 0), 
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which yields the three equations 

-cos a sin 5 z = X 

-sin a sin 5 z = Y 

cos 5 z = Z 


where X, Y, and Z are the components of the vector in inertial coordinates. By convention, 
the angle a is the right ascension of the spin axis, and the angle 5 is the codeclination. The 
declination is 5 = (90° - 5), such that 5 = 90° - 5 substituted in the above equations yields 

-cos a cos 5 z = X; -sin a cos 8 z = Y; and sin 5 z = Z. 

The solution of 5 and a are then given by 

-sin a Y 

= tan a = — 

-cos a. X 


or 



X 
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and 


or 


sin 8 z = Z 



z 


where 


z 2 = X 2 + Y 2 + Z 2 . 


Thus, the spin axis direction can be determined either by finding components in the 
desired reference coordinate system or by solving for the spin axis right ascension and 
declination. 

The remaining constituent of spin-stabilized spacecraft attitude is the rotation around z, 
represented above by 6 , which brings the x-axis into the Zz plane. 


Sensor Angles 

One type of spacecraft sensor measures angles between the spin axis and external physical 
objects or fields, such as the sun or the magnetic field. These sensor angles are related to 
each other in such a way as to make it possible to solve for the spin axis. 

If three different sensors supply the angles, unambiguous solutions are possible. However, 
the use of only two different sensors is preferred, though this produces two-solution 
ambiguities, since one solution can be eliminated by physical arguments, and the savings 
to the spacecraft in space and cost are significant. 


Cone Intersection Interpretation 


The reason that two different sensor angles characterize the spin-axis attitude with a two- 
solution ambiguity, while three sensor angles remove the ambiguity and a single sensor 
gives no useful solutions, is shown by reference to a geometrical argument. (See figure 2 .) 

If the coordinates of a vector? are known in inertial_space, and a sensor angle £ is measured 
between P and the spin axis S, then the direction of S is known to project onto the celes- 
tial sphere at some point on a circle of angular radius £ circumscribing P. 


If two sensor angles, £ and rj, are measured between P and S and between Q and S , respec- 
tively, where? and are known, then the possible orientations of ? are a t the points of 
intersection of the two cones tracing out the circles of radii £ and 17 on the celestial sphere. 
Except for the degenerate case when? and become colinear (within the resolution limits 
of the sensors) and the sensor angles, £ and 17, become equal, at least one intersection of the 
two cones is produced. The one-intersection occurrence is infrequently measured, so that 
a two-solution ambiguity is the general case to be considered in cone-intersection attitude 
determination. 
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Figure 2. Sensor Angle Cones 

One solution can normally be eliminated by considering the past history of the spin axis. 
That is, if the inertial spin axis is a consequence of post-launch procedures, then subsequent 
correct solutions are close to the previous ones, and the unacceptable cone-intersection 
solutions vary rapidly with time. Alternatively, if the time between observations of and 
Q is measured, the ambiguity is resolved since the angular separation of the sensors is 
known. 

After reorientation of the spin axis by ground controlled spacecraft torques, the correct 
solution is again uncertain. Predictions of the effect of the torque should indicate one 
solution as preferable. However, difficulties arise as the intersections draw together, merge, 
and then part. Uncertainty as to the best solution or series of solutions leads to many dif- 
ferent ways of analyzing the sensor angle data and will not be pursued at this time. 

A requirement on the sensor angles due to the motion of the spin axis orientation is that 
they be measured close to the same time— of the order of a spin period. This ensures that 
the actual spin axis will not have moved appreciably between the measurements, such that 
errors between the cone-intersection values and actual spin axis orientations are due pri- 
marily to errors in sensor measurement and telemetry accuracy. 

Finally, if a third sensor is available, then a third sensor angle implies a third cone intersec- 
tion at one of the two previously described cone-intersection loci. This provides a unique 
solution of the spin axis. However, in actual cases, the three cones do not intersect at a 
common locus, but at three close intersections (figure 3). Again, this effect is due to 
errors. 
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Cross Axis Determination 


The cross axis is any axis perpendicular 
to the spin axis; thus, in the above defini- 
tion of spin-stabilized spacecraft attitude, 
the spin axis is the spacecraft z-axis and 
the cross axis is either the x-axis or the 
y-axis. The cross axis is determined by 
specifying the rotation around the z- 
axis which brings the cross axis into a 
known reference plane; for example, 
rotation around z by 6 bringing x into 
the Zz plane. 

To provide known reference-plane/cross-axis relations for the spinning satellite requires 
sensors for observation of external objects or fields, such as the sun or the magnetic field 
of the earth, an internal spacecraft-relative clock or counter of known rate, a means of 
relating the clock to universal time, and the angles from the cross axis to the sensors. 

These requirements serve to relate the cross axis to a direction in space at a particular 
spacecraft-relative time, which corresponds to a particular instant of universal time. 

As an example, consider a spacecraft with a sun sensor (see Appendix B), constant rate 
counters, and associated electronics and telemetry systems. The information for cross axis 
determination is produced through the following procedure: 

1. The angle 6 0 between the cross axis and the sun-sensor/spin-axis plane is con- 
structed according to the spacecraft design. 

2. The spin period T is derived from the number of counts between sequential sightings 
of the sun by the sun sensor. 

3. The number of counts between the transmitting of a particular telemetry word and 
the first sun sighting gives the relative time At t between those events by tj = counts/ 
count rate. 

4. The universal time, t , of the transmission of the particular telemetry word is com- 
puted from the known universal time at the receiving station and the propagation 
time delay from the satellite at known range to the receiving station, Thence, 

5. The universal time of the first sun sighting is 

t = t + At. 

S U 1 

for which time the inertial coordinates of the sun are known. 

6. At time t , the angle 6 in the xy plane between the cross axis and the known 
reference plane is simply the engineered spacecraft angle, 6 = d Q . 
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7. At a somewhat later time, t 8 = t $ + At 2 , the angle is given by 

t 2 

e = e n + — x 360° . 

O lyi 

The final transformation of the cross axis into inertial coordinates depends on the instan- 
taneous position of the sun at time t s in inertial coordinates, and the orientation of the 
spin axis with respect to solar coordinates at time t s , but this initial step in spin-stabilized 
spacecraft attitude determination is made preliminary to and independent of knowledge 
of the spin axis orientation. 

SOME NON-LEAST SQUARES METHODS OF SOLUTION 
The Unique Solution, Three-constraint (PQV) Approach 

If a spacecraft has two cone angle sensors, then at some time t., two vectors,^ and are 
known, and their associated cone angles, j3. and 5 p are measured. The spin axis vector S 
(spacecraft z-axis) is unknown, but is related to the known quantities by the dot products 

A A | 

P. • S - cos j3. = S x + P y S y + P z S z 

i i i 

and 

Q. • S = cos 5. = Q x S x + Q Y S Y + Q z S z 

where P., Q P and S are unit vectors in the XYZ coordinate system. The three spin axis 
components-S x , S y , and S z -are unknown; hence, they cannot be explicitly solved for by 
just these two equations. A third equation in S v , S v , and S, is available, however, by virtue 
of S being a unit vector; that is, 


S x 2 + S Y 2 + S z 2 = 1 . 

This introduction of a quadratic equation leads to a two-solution ambiguity as required for 
the two-cone intersection case. The quadratic approach is algebraically tedious, but a third 
linear equation can be formed which gives a more elegant set of equations for a solution. 
See figure 4. 

A third vector, V, perpendicular to and (f., is defined as the cross product of P. and Q; 

V. = P. X Q. = sin 17 . V. 

1 1 ^1 *11 

where V. is the unit vector. Then, 

„ Pi x Qi 

v. = 

1 sin ri . 
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with components, 



Figure 4. Spin Axis, S, Related to Known Axes P Q., V. 


The angle 17. between P. and Q. is also given by 

P. • Q. = cos 77 . . 

1 ^1 'l 

The relation needed to define the cone angle a. between the new vector V. and the spin 
axis vector S is the dot product 

V. • S = cos a = V X S X + V Y S y + V z S z . 

Thus, to solve the three dot-product equations involving S for the three components S x , 
S Y , and S z requires the evaluation of a. By the use of the law of cosines from spherical 
trigonometry applied to figure 4 and the observation that and equal 90°, the auxiliary 
angles Xj and are determined. From the law of cosines 

cos j3. - cos r) i cos 6 . 

cosxj = — : — 

sin r) i sin 5. 
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and 

A. = 90° - Xj =* cosAj = cos (90° - Xj) = sinXj • 

The desired angle ol now follows, by the law of cosines: 

cos a. = cos 6. cos = sin 6. sin cos A. = sin 5. cos A. = sin 5. sin Xj 
There are now three equations and three unknowns: 

cos 0. = Pv S x + Py Sy + P z S z 

1 1 i i 

cos 5. '= Q x S x + Q y Sy + Q z S z 


COS CL = V x S x + Vy Sy + V z S z 


or in matrix notation 


[Cli = 


P Y. 

1 

P Z. 

l 

q y. 

Qz i 

v Y . 

v z. 


[PQV] j [s]. 


A A 

Unless P is parallel to Q, [PQV] is nonsingular and possesses an inverse such that the solu- 
tion of [S] is 

[S] = [PQV]: 1 [C].. 

The question arises as to what happens to the ambiguity in the initial two-cone-angle solu- 
tion. The choice of the order of P and Q in the cross product makes a difference in the 
sign of Thus, if P X Q = \^, then QXP = - and the resultant [S] matrix components 
from using -\^ in [PQV] are the expected second solutions. From these vector relations 
and knowledge of the spin period, the time between separate sensor sightings, and the angle 
between the two sensors, an unambiguous method of selecting the correct solution is possi- 
ble. This method and appropriate tests are described in Appendix C. 

If the time between sensor sightings is not available, as is the case with magnetometers, 
many sets of equations can be formed corresponding to different conditions at different 
times t., with each set of equations providing two solutions as outlined above. From each 
pair of solutions, one is closest to previous solutions or a priori expectations of the spin 
axis vector, and is chosen for the set of all best solutions. These best solutions will be 
reasonably close in value, but will differ even for solutions of a stationary spin axis vector. 
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The errors which propagate to create the solution discrepancies are due (1) to engineering 
limitations (for example, half-degree accuracy of the sun sensor, three quarters of a degree 
accuracy for an earth sensor, and angle uncertainties related to the digital counter); (2) to 
field model uncertainties; (3) to spacecraft, sun, and other position vector uncertainties; 
and (4) to random errors at each stage of the procedure. The systematic errors can be 
lessened by various long-period, observation-deduced adjustments, but there are always 
residual random deviations. 

Three-constraint Solution Applied to IMP-6 

IMP-6 is the sixth in the series of Interplanetary Monitoring Platforms. It has a highly 
elliptical orbit about the earth (perigee of approximately 12,000 km and apogee of 
approximately 210,000 km), and is primarily engaged in investigation of the earth’s 
magnetic field. For attitude determination, it possesses an earth horizon sensor and a sun 
sensor. 

The PQV method for IMP-6 attitude determination makes use of the satellite’s earth hori- 
zon and sun sensors to derive cone angles. The derivations of the earth nadir cone angle 8 
depend on the earth-spacecraft-sun relations, with two cases favorable for cone angles. 

Case 1 is involved when the earth sensor detects two earth horizons in a single sweep across 
the visible earth. Case 2 applies when one horizon and one terminator are detected. 

The sensors operate generally as described in Appendix B. Specifically, a 9-bit sun sensor 
gives 0.5° accuracy in the sun angle, )3; a 1.5° horizon sensor field of view gives a 0.75° ap- 
I parent increase in the earth radius; and a 1600 counts-per-second oscillator drives the count- 
ers and thus provides resolutions of 0.000625 s between counts. The counters are a sun 
clock measuring the time between a specific bit of the telemetry format and the next sun 
pulse, a spin-period counter measuring the time between successive sun pulses, a sun-earth 
horizon counter measuring the time between a sun pulse and the next initial sensing of the 
earth, and an earth-width counter, which measures the scan time between initial earth hori- 
zon detection and the trailing earth horizon. The mounting angle 7 (figure 5) between spin 
axis and earth sensor is 87°, off by 3° from the desired 90°. 

The presence of a terminator, rather than an earth horizon, as either the initial or trailing 
• counter trigger is determined analytically. From the sun-earth horizon counter, El (Earth 
IN), and the spin period counter, SPC, the angle v between the spin axis/sun plane and the 
spin axis/earth horizon plane is found by 

EIT = El counts X 0.000625 seconds/count 

SP = SPC counts X 0.000625 seconds/count 

v EIT 0 EIT 

— = or >- = 360 X 

360 SP SP 
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Figure 5. Vectors and Related Angles 


where EIT 

= Earth In Time, 

El 

= Earth In, 

SP 

= Spin Period, 

SPC 

= Spin Period Counter. 


If 0 < v < 180°, initial earth detection is at the earth horizon, because of the sunlight-on- 
earth illumination geometry. The angle 2 p subtended by the earth at the spacecraft distance 
is derived from sin p = R e /(R e + h), where R e is the earth radius (see figure 6a) and h is the 
spacecraft altitude. 

The angle 0 between the spacecraft position vector and the sun vector is defined by 
P • (-Q) = cos 0 or by 


< p = 180° - (u + p) for 0 < v + p < 180° 

<p = (v + p) - 180° for 180° <v + p < 360° . 
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Figure 6. Relation of Spacecraft to Sunlit and Shadowed Earth 


Then, for 0 < p, as for the spacecraft in position 1 in figure 6b, the spacecraft views a 
fully sunlit earth, and the earth sensor senses two earth horizons. But for 0 > p and v > 0, 
with the spacecraft in position 2, the spacecraft senses an earth horizon and a terminator. 
When v = 0, with the spacecraft in position 3, then either the earth horizon sensor does not 
scan the earth or the earth viewed by the sensor is entirely shadowed. 

With these criteria established, the earth nadir angle, 5, can be derived according to one of 
the two cases mentioned above. 


Case 1— Nadir Angle Derivation for Scan of Two Horizons— In the event that 0 <p (figure 
6b, spacecraft at position 1), a scan from horizon to horizon by the earth horizon sensor 
(figure 7) provides a measure of an arc, p, across the earth. Half of the angle p is the angle 
p/2 between the spin axis-earth horizon plane and the spin axis-earth center plane. By 


the law of cosines, 
and by 


cosp = cos 7 cos 5 + sin 7 sin 6 cos — , 


cos 2 5 + sin 2 5 = 1 


the equation can be solved. Another linear relation exists, however, and can be used rather 
than the quadratic: 

fl 

cos 17 = cos |3 cos 5 + sin |3 sin 6 cos (V + — ■ ) . 
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Figure 7. Case 1, Nadir Angle Derivation, Scan of Two Horizons 

The cos 77 and cos p equations can be viewed as two equations in the two unknowns, cos 5 
and sin 5 , and solved accordingly: 


such that 


cos 5 


cos p - sin 7 sin 5 cos— 

2 

cos 7 


cos 77 cos 7 - cos 7 sin 0 sin 5 cos 


KH 


mV 

cos p - sin 7 sin 5 cos — 1 cos j 3 sin 5 
cos p cos j 3 - cos 77 cos 7 sin 5 


[■ 


cos 0 sin 7 cos — - cos y sin 0 cos 
2 


(">)]■ 


cos 0 cos p - cos 77 COS 7 


cos 0 sin 7 cos— - cos 7 sin 0 cos 


K> 


All angles on the right-hand side of the equation are determined from spacecraft measure- 
ments ((3, and p/2), from orbit determination (p and 77), and from spacecraft engineering, 
(the mounting angle of the horizon sensor 7), allowing 5 to be computed. 
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Case 2— Nadir Angle Derivation for Scan of One Horizon and the Earth Terminator— In the 
event that <f> > p and v > 0 , no useful information is acquired from the earth horizon sensor 
scan from an earth horizon to the earth terminator (figure 8). The angle g merely provides 
an angle that is less than the full arc across the earth. Elaborate geometrical and trigono- 
metric analysis could supply a proportionate angle equivalent to g/2, but a simpler approach 
is available. 


A 


S 



Figure 8. Case 2 , Nadir Angle Derivation, Scan of One Horizon and Earth Terminator (7<6) 

The great circle arc X, between the sun position vector P and the earth horizon vector ft, is 
found by use of the law of cosines, 

cos X = cos 0 cos 7 + sin 0 sin 7 cos v . 

The interior angle £ results from an application of the law of sines, 

sin v 

sin £ = sin 7 . 

sin X 
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The last auxiliary angle, the interior angle e, is also given by the law of cosines: 


Then, 

for 7 < 5 and 


cose = 


cos p - cos X cos 17 
sin X sin 77 


cos 5 = cos j 3 cos 77 + sin j 3 sin r\ cos (£ 


cos 5 = cos 0 cos 77 + sin (3 sin 77 cos (£ 


+ e) 

- e) 


for 7 > 5 (figure 9). But since 5 is not known beforehand in relation to y, a two-value 
ambiguity in 5 results. 



Figure 9. Case 2, Nadir Angle Derivation, Scan of 
One Horizon and Earth T erminator (y > 5 ) 


Various other angles can be developed 
and applied as criteria to the selection 
of allowable data for the solution of 5. 
(For this more extensive and complete 
treatment of the PQV method see 
reference 3.) 

One other augmentation of the above 
equations is necessary for completeness. 
If the terminator is sensed first by the 
earth scanner, and an earth horizon 
triggers the trailing pulse, then in each 
of the equations for case 2 where v is 
present, the substitution v + n is made. 


A relation mentioned above in another context also provides the specifications for the order 
of the horizon and terminator. If 0 < v + p < 1 80°, an earth horizon is sensed prior to the 
terminator. If 180° < v + p < 360°, the earth terminator preceeds the horizon. For v + p 
angles near 180° (dependent on the radial distance of the spacecraft from the earth), the 
spacecraft views a fully sunlit earth. 


By these means the required cone angles j3 and 5 are acquired for the solution of the IMP-6 
spin axis orientation. The solution itself follows the algorithm outlined above with the 
additional ambiguity in 8 in the horizon/terminator case leading to a four-way ambiguity 
in the spin-axis solution (two 5’s for each cone intersection). 

An adjustment of the 6 and ju angles must be made in practice in order to eliminate the 
systematic error reducing the counts caused by the 1 .5° field of view of the telescope. 

This adjustment is discussed in the last paragraphs of “Spherical Coordinate Solutions Applied 
to IMP-6,” below. 
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The Unique Solution, Spherical Coordinates (Right Ascension and Declination) Approach 

A more general treatment of the spin-stabilized spacecraft attitude problem focuses on the 
long-term, inertially-constant axis of rotation. This axis, the net angular momentum vector 
£, is not, in general, the instantaneous spin axis S . Except for the special case where the 
principal moment of inertia I z of the z-axis is greater than the transverse principal moments 
of inertia I and I , the angle between the two moments at which the spin axis precesses 
about the momentum vector tends to zero (unless the symmetry axes and principal moment 
axes are misaligned). This is the nominal condition of most spin-stabilized spacecraft, and 
the following discussion makes reference to the momentum vector, rather than to the spin 
axis, in accordance with existing literature* on the method (reference 4). 

Figure 10 is a representation of the celestial sphere centered at the spacecraft center of 
mass with all vectors emanating from that point. The components of the vectors are in 
spherical coordinates with respect to the equatorial, inertial reference system defined by 
Z = celestial north, X = first point of Aries, and Y = orthogonal axis making a right-handed 


z 



Figure 10. Right Ascension and Declination Components of Earth, Sun, and Momentum Vectors 


♦Pyle, E. J. and J. S. Albus, IMP F & G Aspect Determination System , GSFC X-71 1-67-30 (TM X-63526), January 1967. 
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system. The spherical components of a unit vector V are the right ascension RA v and the 
declination D y . The vectors shown are the sun position, earth nadir, and angular momentum 
unit vectors, 2, N, and L, respectively. The arcs of the great circles between the vectors 
are 0 between 2 and L, rj between 2 and N, and 5 between L and N. 

As discussed in the section, “The Unique Solution, Three-Constraint (PQV) Approach,” 
the angles )3 and 8 are determined from spacecraft observations, 2 is the sun position at the 
time of the observations, N is the negative vector of the spacecraft orbital position, and t? 
is the arccosine of the dot product of 2 and N. From these known quantities, the angular 
momentum components, RA L and D L , can be derived. Explicitly, use is made of the 
spherical coordinates, RA 2 , D 2 , RA N , and D N , the angles 5, )3, and ij, various auxiliary 
interior angles, and the spherical trigonometric laws of sines and cosines in the following 
equations to obtain the required relations. 

Let the auxiliary interior angle 0 be formed by the ccw rotation of the 5 arc into the rj arc 
around N. The angle is defined by the law of cosines relation 

cos 0 - cos t? cos 6 

cos 0 = 

sin t? sin 5 

Similarly, the auxiliary interior angle between (90° - D N ) and r? is formed by 

cos (90° - D £ ) - cos (90° - D n )cosi? sin D 2 - sinD N cost? 

cos 9 = - = . 

sin (90 -D N )sinr? cos D n sin t? 

The difference angle (0 - 0) is used to define the D L relations, 

cos (90° - D L ) = cos (90° - D n )cos 5 + sin (90° - D N ) sin 5 cos (0 - 0) 

and 

sin D l = sin D N cos 6 + cos D N sin 6 cos (9 - 0) (-90° < D L < 90°). 

NOTE: If t? is between 5 and (90° - D N ), then 0 becomes a negative angle by 

the defining convention. Sign considerations must be implemented 
on the D l equation for full generality. This is done by using the 
additional defining equations of 0 and 6 to provide the quadrant of 
each angle. 

A second 0 relation is 

sin £2 

sin 0 = sin 0 

sin r? 

where SI is the interior angle about L from the 2 L plane to the NL plane. This angle is 
determined by sensor pulse-times, and will be discussed in the next section, “Spherical 
Coordinate Solutions Applied to IMP-6.” The second d relation is 



sin (RA n - RA S ) sin(RA N - RA £ ) 

sin 0 = sin (90° - D„) = cos D_ 

^ sin t? ^ sin r? 


Two relations are also required to define RA L and its quadrant 

cos 6 - cos (90° - D n )cos(90° - D L ) cos 6 - sin D N sin D L 


cos(RA n - RA l ) = 


sin (90° - D n ) sin (90° - D L ) 


cos D n cos D l 


and 


sin (RA n - RA l ) = sin 5 
Since RA N is known, 


sin (0 - <{> ) 


sin 5 


sin (90 - D l ) 
ra l = RA n - (RA n - RA l ) . 


sin (0 - <t>) 
cos D. 


The parameters given in the section, “The Unique Solution, Three-Constraint (PQV) Ap- 
proach,” are equivalent to those of this section by 



Figure 11. Right Ascension, Declina- 
tion and XYZ Components of L 


P o i 
Q o N 
S ♦ L 

and the spherical coordinate components of L, as 
shown in figure 1 1, are related to the Cartesian coor- 
dinate components of S by 

X = cos D. sin RA. 

s L L 

Y = cos D t cos RA t 

S L L 

Z e = sin D r . 

S L 


Spherical Coordinate Solutions Applied to IMP-6 

The application of the RA and D method to the IMP-6 satellite obviously entails the use of 
the same telemetry information produced by the same engineering design described above 
in the section, “Three-constraint Solution Applied to IMP-6.” However, different algorithms 
for the determination of the angle 5 between the L and N vectors were incorporated in the 
actual RA and D method attitude determination program. Accordingly, these algorithms 
and the pertinent spacecraft observations are briefly examined below. 



First, consider case 1 for which two earth horizons are scanned. In this case, the angle /x 
between the horizons is determined from the earth width counter, and the angle v between 
the 2L plane and the LH r plane is determined from the sun-earth counter. The angle 
/x/2 between the LH l plane and the LN plane provides the final angle necessary to form 
the interior angle £2 of the previous section, £2 = v + /i/2 (figure 12). 



Figure 12. Case 1, Two Earth Horizons Scanned 


By the law of cosines. 


cos r? = cos /? cos 6 + sin /? sin 5 cos £2 


and by 


sin 5 



cos 7 } 


cos /? cos 5 = 


sin /? cos £2 


y/l~ - cos 2 5 . 


Squaring both sides and collecting terms of order cos 5, 

(cos 2 j3 + sin 2 /? cos 2 £2) cos 2 6 - (2cos r? cos j8) cos 5 + (cos 2 r? - sin 2 /? cos 2 £2) = 0. 
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By the quadratic formula, 


cos 6 = 


2 cost? cos /3 ± v^cos^rjcos^j3-4(cos 2 0 + sin 2 /} cos 2 J2) (cos 2 tj - sin 2 /} cos 2 12) 
2 (cos 2 /} + sin 2 /3 cos 2 £2) 


which simplifies to 


2 cos?] cos /} ± sin /3 cos Cl \/cos 2 /3 + sin 2 /} cos 2 i2 - cos 2 i? 

cos 6 = 

cos 2 /} + sin 2 /} cos 2 Cl 


This two-solution ambiguity in 6 is resolved by consideration of whether the angle 5 should 
be less than or greater than the mounting angle 7 of the sensor with respect to the spin axis. 

Next, consider case 2 for which one earth horizon and the earth terminator are scanned 
(figures 13, 14). The items of significance in this case are the angles jS, 7?, 7, v, and p, all 
defined previously, and the auxiliary angles, X, B, D, and E, derived from the known. The 
terminator /scan intersection T provides no information applicable to this derivation. 



Figure 13. Case 2, Earth Horizon and Earth Terminator Scanned 
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Figure 14. Earth Horizon and Earth Terminator 
Scanned (Detail of Figure 13). 


From the law of cosines, 
cos X = cos (3 cos y + sin P sin y cos u 

which holds for 0° < v < 180°. 

For v > 180°, the terminator is 
scanned first and the horizon, 
second, requiring v to be replaced 
by ju + v in the equation. If v = 1 80°, 
then X = /3 + 7. For this case, v is 
the angle determined by the space- 
craft counter which is triggered by 
the sun and shut off by an earth 
horizon/terminator scan, and p is 
derived from the counter which is 
triggered by the horizon/termina- 
tor and shut off by the terminator/ 
horizon. 


The interior angles B, D, and E follow from the application of the laws of sines and cosines. 

cos /3 - cos 7 cos X 


cos B 


sin 7 sin X 


and 


sin B = sin j 3 


sin v 


sin X 


For the interior angle E, two equations are not available to define the angle and its quadrant; 
so, a two-possible-values cosine equation produces an ambiguity in E, 

cos 77 - cos X cos p 

cos E = , 

sin X sin p 

where E can be in either of two quadrants. As a result, there are two possible values of D 
determined by 

D = 360 ° - B - E 

and the defining equation for 5 in this case propagates a two-value ambiguity, 

cos 6 = cos p cos 7 + sin p sin 7 cos D, 
where 5 in either of two quadrants satisfies. 

The desired orientation for IMP-6 at the time of launch was for RA L = 90°, D L =-66?5, 
and a spin period of 1 1 seconds. After a series of maneuvers, from launch on March 13, 

1971 , through July 15, 1971 , the satellite was permitted to remain in a spin-stabilized, 
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spin-axis orientation of about RA L = 95°, D L = -67?5, and a spin period of 1 2.8 seconds. 
In the course of determining the definitive attitude, an apparent coning of the spin axis 
about a separated momentum vector was not a necessary hypothesis if an assumption was 
made that the earth horizon sensor mounting angle 7 was 87° rather than the desired, 
engineered 90°. 

Another value adjusted to bring all solutions over an orbit into closer relative agreement 
was the angular earth radius p. The reason for this adjustment is the triggering of the earth- 
horizon scan pulse before the center of the sensor reaches the horizon (figure 1 5). The 
approximate adjustment value is half the scan width r, that is, the radius of the field of 
view of the sensor telescope. Thus, when a certain small percentage of the sensor field of 
view is covered by the sunlit earth, the first horizon pulse or the trailing pulse is initiated. 
This provides fewer counts to the sun-earth horizon counter and more counts to the earth 
width counter, in direct proportion to the adjustment in p. The best value for the adjust- 
ment shifts according to the spacecraft distance to the earth, because the relative size of 
the earth with respect to the scan width varies, and the greater the curvature, the greater 
the proportion of the sensor which must cross the horizon to trigger the pulse (figure 1 6). 



Figure 15. Earth Horizon Sensor at Times of Scan Pulses 

For IMP- 6 , with a nominal scan width 1 ? 5 , the adjustment value near perigee has been 
0?8 and the value near apogee, 0?6. The chosen adjustment value affects only those 
equations with p related to counter-derived values. 
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(b) Spacecraft at Apogee 


Figure 1 6. Different Adjustments of Angular Earth Radius p with 
Constant Percentage of Sensor Covered 


LEAST SQUARES METHODS OF SOLUTION 


The Iterative Differential Correction Technique 

In a typical spin-stabilized spacecraft attitude determination system an abundance of cone 
intersection pairs are produced over a relatively short time by the spacecraft sensors. To 
make use of all of them, the method of least squares is applied to solve the over-determined 
set of equations. See Appendix D for a brief description of the method of least squares. 

If the equations involve spherical coordinates, and are to be solved for the angular attitude 
components RA L and D L , then the equations are not linear and a closed form solution in 
RA l and D L is not practical. Instead, a series of small corrections to successive, iterated 
values of RA L and D L are produced, until the size of corrections are less than preselected 
small amounts. 


Each set of solution correction values is produced by a least squares treatment. The pro- 
cedure is to expand the equations in a Taylor’s series in terms of the required coordinates, 
truncate the series after the linear terms, solve the truncated equations by least squares for 
the correction values, apply the corrections to the last values of the components, and 
repeat the cycle until the newest corrections are less than the preselected bounds. 

The form of the equations prior to expansion in terms of RA L and D L is, for n equations, 

U. • L = cos 6 i = U x cos D l cos RA l + U y cos D L + RA L + U z sin D L 


where U. 

i 


L 

A; 

u x ..u Y .,u z 


the known unit vectors (sun, earth nadir, or magnetic field), 
the spin axis (or angular momentum) unit vector, 
the cone angles (p, 8, and so on), 

the inertial coordinate components of U., i = 1 , . . . n, and 
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D l and RA l = the inertial, spherical coordinate components of L, which are 
to be solved for. 

The functional equation 

Fj (RA l , D l ) = cos 0. - U x cos D l cos RA l - U y cos D L sin RA L - U z sin D L = 0 


is expanded in a Taylor’s series about initial coordinates RA° , D° by 


/ 3F. \ 

F i (RA L> D L ) = F i (RA° L> D° L ) + 


ara l + 


' 3F. 
3D 


ad l + 


L/ 0 


'0 2 F. 


3RA 2 


(ARA) 2 + . . . 


and truncated after the (3F j /9D L ) 0 AD L term. For notational convenience let 


then 



3F. 


3RA: \ 3D 



3F. 


3D; 


3F. 


3F. 


Fj (RA l , D l ) - Fj (RA° , D° ) + ARA l + A D L 


3ra: 


3D, 


where 


Fj (RA° , D° ) = cos 0j - U x cos D° cos RA° - U y cos D° sin RA° - U z sin D° = p° 


and the p° are the residuals of the initial estimate. In matrix notation and rearranging 
terms, noting that F. (RA L , D L ) = 0, 


Pi 


~3Fj/3RA° 3Fj/3D°~ 


P 2 


3F 2 /3RA° 3F 2 /3D° 


ARAj 

. 

= 

• • 

• * 


- AD l. 

_Pn_ 


• • 

_3F n /3RA“ 3F n /3D°_ 


(nXl) 

(nX 2) 

(2 X 1) 


= R = PA. 


If an n-by-n diagonal weighting matrix K' 1 can be constructed from known error magnitudes 
in the sensor’s cone angle measurements, such that the relative accuracies of each measure- 
ment can be reflected in the solution, then 


K 1 R = K 1 PA 
(n X 1) (n X 1) 
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Specifically, 


1 

o 

• 

• 

• 

o 

1 


a 2 (cos 0^0 * • • 0 

0 k 2 • *0 

• o • • 


0 o 2 (cos0 2 ) • * • 0 

* • • • 

1 

o 

o 

• 

• • 

j* ’ ‘ 

1 


« • • • 

• • • • 

0 0 • o 2 ( cos 0 n ) 


where o 2 (cos 0 ; ) = +sin 2 0.o 2 in analogy with A(cos 0) = -sin 0 A (0). 

u i 

The solution for A now follows the usual matrix least squares procedure. First, premultiply 
both sides of the equation by P T . 

P T K*R = P T K 1 P A ^ 

(2X1) (2X2) (2X1) 

If the inverse of (P T K' 1 P) is nonsingular, 

(p T K -l pyl pT^l R = (pT K -lp)-l (pT K -l p) A = A 

which is the initial correction matrix with elements ARA° , AD° . The correction elements 
are combined with the initial estimates RA^ , D° to form the next estimates in the iteration 
process. Before returning to the Taylor’s expansion, ARA° and AD° are compared with 
the preselected bounds and e 2 . If 

ARA° < e 1 and AD° < e 2 , 


then the solution of the spin axis orientation is accepted as 


RA° + ARA° and D° + AD° . 

If, as is usual, one or both of the bounds is exceeded, the algorithm repeats. 

The inverse of the normal equations (P T K' 1 P)' 1 for each iteration provides another desirable 
set of solutions; that is, (P T K ‘ P)' 1 is the 2-by-2 variance-covariance matrix of error in the 
solution correction values. 


(p T K -l pyl 


o 2 RA l oRA l oD l 

aRA L aD L ° D L 2 


The matrix is formed at each iteration, and when the correction values are found to be 
acceptable, the standard deviations <*RA L and aD L are available. 
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Table 1 

Conversion of Digital Solar Attitude Detector 
Gray Codes to Sun Angle* 


Code 
Bit No. 
7654321 

Angle X 
(deg) 

Code 
Bit No. 
7654321 

Angle X 
(deg) 

Code 
Bit No. 
7654321 

Angle X 
(deg) 

1110101 = 

0.5 

1110100 = 

1.5 

1111100 = 

2.5 

1111101 = 

3.5 

1111111 = 

4.5 

1111110 = 

5.5 

1111010 = 

6.5 

1111011 = 

7.5 

1111001 = 

8.5 

1111000 = 

9.5 

1101000 = 

10.5 

1101001 = 

11.5 

1101011 = 

12.5 

1101010 - 

13.5 

1101110 = 

14.5 

1101111 = 

15.5 

1101101 = 

16.5 

1101100 = 

17.5 

1100100 = 

18.5 

1100101 = 

19.5 

1100111 = 

20.5 

1100110 = 

21.5 

1100010 = 

22.5 

1100011 = 

23.5 

1100001 = 

24.5 

1100000 = 

25.5 

0100000 = 

26.5 

0100001 = 

27.5 

0100011 = 

28.5 

0100010 = 

29.5 

0100110 = 

30.5 

0100111 = 

31.5 

0100101 = 

32.5 

0100100 = 

33.5 

0101100 = 

34.5 

0101101 = 

35.5 

0101111 = 

36.5 

0101110 = 

37.5 

0101010 = 

38.5 

0101011 = 

39.5 

0101001 = 

40.5 

0101000 = 

41.5 

0111000 = 

42.5 

0111001 = 

43.5 

0111011 = 

44.5 

0111010 = 

45.5 

0111110 = 

46.5 

0111111 = 

47.5 

0111101 = 

48.5 

0111100 = 

49.5 

0110100 = 

50.5 

0110101 = 

51.5 

0110111 = 

52.5 

0110110 = 

53.5 

0110010 = 

54.5 

0110011 = 

55.5 

0110001 = 

56.5 

0110000 = 

57.5 

0010000 = 

58.5 

0010001 = 

59.5 

0010011 = 

60.5 

0010010 = 

61.5 

0010110 = 

62.5 

0010111 = 

63.5 

0010101 = 

64.5 

0010100 = 

65.5 

0011100 = 

66.5 

0011101 = 

67.5 

0011111 = 

68.5 

0011110 = 

69.5 

0011010 = 

70.5 

0011011 = 

71.5 

0011001 = 

72.5 

0011000 = 

73.5 

0001000 = 

74.5 

0001001 = 

75.5 

0001011 = 

76.5 

0001010 = 

77.5 

0001110 = 

78.5 

0001111 = 

79.5 

0001101 = 

80.5 

0001100 = 

81.5 

0000100 = 

82.5 

0000101 = 

83.5 

0000111 = 

84.5 

0000110 = 

85.5 

0000010 = 

86.5 

0000011 = 

87.5 

0000001 = 

88.5 

0000000 = 

89.5 


*If bit 8 (not listed above) is a 1-bit, the sun angle = angle X. 

If bit 8 is a 0-bit, sun angle = 1 80° - angle X. 

Bit outputs other than those shown result in a flag in the angle value. 
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where H x , H y , and Hj, are the measured components. 

The field components are read out in analog form in millivolt units from -250 to +250 milli- 
volts, which correspond to physical units of -27.853 to +27.853 amperes/meter (-350 to 
+350 millioersteds). For telemetry back to earth, the millivolt readings are digitized as 
8-bit words. But 8 bits can represent just 256 readings, so that digitization introduces a 
±1 millivolt error in each measurement. The formula for converting from digital counts 
back to millivolts is as follows: 

If counts < 127, millivolts = 2 X (counts - 127). 

If counts > 127, millivolts = 2 X (counts - 128). 

This shows that two different counts for zero millivolts is possible-127 or 128. 

The final conversion from millivolts to millioersteds is straightforward, since the two values 
are proportional. Thus 

700 

Millioersteds = (millivolts + 250) X 350 

500 

where (millivolts + 250) = converted millivolts minus lower millivolt limit, 

700/500 = the ratio of the range in millioersteds to the range in millivolts, 

and 

-350 = the addition of the lower millioersted limit. 

This formula is a linear interpolation technique. 

Besides the digitization error of ±1 millioersted, a small error in the magnitude of the mag- 
netic field is introduced by the sampling of the three magnetometers over three consecutive 
half frames of 0.768 second once each major frame of 49.152 seconds, rather than at the 
same time. The total effect of all magnetometer measurement errors is estimated to be 
±280 milliam peres/meter (±3.5 millioersteds) for the system under operating conditions.* 

A chargeable, trim-magnet system, of three chargeable magnets parallel to the spacecraft 
axes, is included to provide a means of compensating for residual dipole fields from the 
z-coil if demagnetization is not fully effective. These three trim readouts are in analog 
voltages corresponding to pole-cm units. The pole-cm units are, in turn, proportional to 
millioersted units, and are then applied as residual bias corrections to the respective mag- 
netometer values. 

Additional analyses of magnetometer data were performed to determine and eliminate other 
types of systematic error in the magnetic field/magnetometer cone intersection angles. The 
sources may have resulted from launch effects or other post-launch conditions, hence 


*Meyers, G. F., et al., SAS-1 Attitude Data Analyses, GSFC X-542-71-363, 1971. 
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the systematic errors had to be traced from in-flight information. The principal possible 
sources of systematic errors were sensor mounting misalignments, magnetometer non- 
orthogonality, and inate sensor biases. Two different kinds of mounting misalignment 
were considered: (1) If the spacecraft geometrical symmetry axis and spin axis produced by 
the spin rotor do not coincide, but are some small angle apart, then the magnetometer 
triad is misaligned to the spin axis, even when the triad is mounted correctly to the sym- 
metry axis. (2) If the symmetry and spin axes coincide, but the tri-axis magnetometer is 
mounted with the magnetometers not parallel to their respective spacecraft axes, then the 
triad is misaligned to the spin axis. 

If the individual magnetometers are mounted improperly with respect to each other, such 
that they are not mutually perpendicular, then systematic errors due to nonorthogonality 
arise. Finally, if after the chargeable trim magnets and demagnetization of residual coil 
fields are considered, and magnetic field components do not match the spacecraft magne- 
tometer measurements within a random range, then a remnant, residual spacecraft magne- 
tism and sensor bias is present. Such a residual bias can be present separately for each mag- 
netometer. 

The chief difficulty in determining each of the above factors is that they are interdependent. 
As such, for instance, a misalignment error can easily be interpreted as a contribution to 
the residual bias. Cross-check and multiparameter programs were developed to circumvent 
the interdependence problem. The best determination of each of the factors is given in 
tables 2, 3, and 4. Note that the misalignment factors are not separated according to the 
two kinds of spin axis-to-magnetometer misalignment. A further study of short term drift, 
caused by either kind of coning of the tri-axis magnetometer about the spin axis, indicated 
a possible coning amplitude of about 0?6 with a period of the order of a spin period. How- 
ever, the two kinds of misalignment produced the same effects in simulated data, and noise 
effectively obscured any oscillation attributable to a specific misalignment. 

In case of telemetry error or other kind of data loss from one of the magnetometer com- 
ponents, the magnetic field cone angle n can still be determined from the remaining com- 
ponents, and the geomagnetic field magnitude can be taken from special orbit data tapes 
or computed from a programmed magnetic field model. If H is the magnetic field magni- 
tude in milligauss and M x , M y , and M z are the field components derived from the magne- 
tometer measurements, then cone angle ju is given by either of the following equations: 


cos jU 


x/h 2 - M x 2 - M y 2 

for M not known. 

H z 


cos n 



for M x and/or M y not known. 
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Table 2 

Angles Between Magnetometers 


Between 

Angle (Degrees) 

x and y* 

90.775 ± 0.090 

x and z 

90.19 ±0.38 

y and z 

89.55 ±0.47 

* Angle not corrected for 0?37 rotation of spacecraft between x and y readings. 

Table 3 

Alignment of Magnetometer to Spin Axis 

Magnetometer 

Angle from Magnetometer 
to Spin Axis (Degrees) 

X 

90.06 ±0.41 

y 

89.75 ±0.49 

Z 

0.307 ± 0.088 


Table 4 

Magnetometer Residual Biases 


Magnetometer 

Bias 

Milligauss 

Nanotesla 

X 

-0.23 ± 0.57 

- 23 ± 57 

y 

3.21 ± 0.30 

321 ± 30 

z 

6.23 ± 0.94 

623 ± 94 


If all components are known, the relation is 

M z 

COS J u = — 

Vm x 2 + M y 2 + M z 2 
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The Closed Form Least Squares Method 

Another least squares treatment of the spin-axis orientation problem, but one which does 
not require iterations with small corrections to an estimated solution, is the closed-form 
least squares (CFLS) method. From the initial set of cone angle equations, one pass through 
the solution algorithm derives the results. However, an exception to the one-pass condition 
exists if some of the cone angle data is bad enough to be discarded, and the algorithm is 
repeated with the reduced data set. Various acceptance criteria could be selected to desig- 
nate bad data, but the most common is to choose a specific number of standard deviations 
from the solution parameters as the acceptance limit. 

The solution components are the spin-axis inertial Cartesian coordinates, X s , Y s , and Z s , 
which are involved through the dot-product/ cosine-of-the-cone-angle relation 

cosfl; = U. • S = U x X s + U y Y s + U z Z s (i = l,2,...,n) 

where 0. is the cone angle between the i th physical, unit vector U. and the spin axis S. The 
n equations of this form can be expressed in matrix notation as 


(nX 3) 


cos 0 j 



\ 


cos 0 2 


% 

U v, 


• 


• 

• 

• 

• - 


• 

• 

• 

♦ 


• 

• 

• 

COS d n 


U x 

n 

u Y 

n 

U Z 

n 


(3X1) 


X 


which, for convenience, is represented by 

C = U S (condition equations). 

Using the matrix least squares techniques of Appendix D, premultiply by the transpose of U 

U T C = U T U S (normal equations) 

(3x7) (3x7) (3X1) 

where U T U is a 3-by-3 square matrix. If U T U is nonsingular, then it has an inverse and 

[U T U]' 1 U T C = S 

where the elements of S are the components X s , Y s , and Z s . If required, the spherical coor- 
dinates are derived as in the section, 4t Spin-stabilized Spacecraft Attitude, and are illus- 
trated in figure 1 1 . 
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The addition of a weighting matrix is desirable if a priori relative sensor precision is 
unknown. But, the superiority of one sensor over another is reflected in the standard 
deviations, a ; , of each cos d. measurement. The weighting matrix K 1 is therefore given 
by the inverse of the diagonal variance matrix K of the measurements, where o. 2 = variance 
and sin d.o a - a,. 

1 tt . I 

1 

o t 2 0 • • • 0 

0 o 2 2 • • • 0 


K = 


0 0 • 


2 


The initial equations then take on the form 

K ' 1 C = K ' 1 U S (weighted conditional equations), 

U t K -1 C = [U t K‘* U] S (weighted normal equations). 


The 3-by-3 inverse of [U T K 1 U] in the normal equations is the variance-covariance error 
matrix for the solution components from which the standard deviations o x , o Y , and o z 
are produced. 


[l^K ’ 1 U ] 1 


a X a X a Y °X °Z 

CT Y a X CT Y 2 °Y a Z 
°Z °X °Z °Y °Z 


Closed Form Least Squares Solution Applied to IMP-6 

The IMP-6 spacecraft attitude system and a pair of exact, three-constraint approaches to 
spin-axis determination have been described above. Application of the CFLS method paral- 
lels that of the PQV approach for the computation and selection of the earth nadir/spin 
axis angle 5. Because the CFLS method is applied to overdetermined sets of data, additional 
equations such as for the V vector or the unit spin axis vector S constraint are not required 
to produce a solution. However, evaluation of the V vector components and the angle a 
between V and S (as done for the PQV approach) provides a third set of equations which 
are combined with the others for a better least squares solution. The advantage of adding 
the V-terms is to constrain the solution in such a way that errors in the P and Q terms can- 
not combine to bias it in the direction of the V’s. 
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If measured cone angles from the spin axis to the sun and to the earth nadir varied widely 
between 0° and 180°, constructed vector constraints would not be necessary; but, over the 
period of time for which the spin axis can be considered stationary, the sun position changes 
little, and the earth nadir direction and cone angle 5 are useful only when the earth is 
scanned by the earth horizon sensor. Thus, different cone angles in the V-terms serve to 
improve the least squares solution of S. 

The condition equations for the CFLS solution of the cone-intersection attitude problem 
in general form is therefore 

cos 9 ^ = U x X s + U y Y s + U z Z s 0 = l,2,...,3n) 
where the index j extends to 3n to represent the three specific condition equation subsets: 


cos 4 = P x X s + P y Y s + P z Z s 

0 = 

1,2,.. 

• , n), 

cos Sj = Q x _X s + Q Y .A + QzA 

0 = 

1,2,.. 

• , n). 

COSO. = V x X s + V y Y s + V z Z s 

0 = 

1,2,.. 

. . , n), 


where the index i is related to a specific telemetry readout. 

ANALYSIS AND EVALUATION 
SAS-1 Attitude Solutions 

As described above, SAS-1 spin-axis solutions are solved by the cone intersection IDCT 
approach from magnetometer and sun sensor data. Twenty-eight early orbits of SAS-1 were 
analyzed by the Attitude Determination and Control Section, GSFC.* These orbits were 
selected because they occurred before the spacecraft tape recorder stopped functioning. 

Thus an abundance of good playback data from an entire orbit was available. From these 
analyses several additional parameters were obtained, such as the orthogonality charac- 
teristics of the magnetometer axes, the misalignment of the magnetometer axes to the 
spacecraft axes, and uncorrected residual magnetic field biases. Nonorthogonality and 
misalignment , adjustments were shown to be slight, but the residual bias corrections im- 
proved the solutions appreciably. 

For this report, raw data from 1 1 of the 28 SAS-1 analyzed orbits were acquired and were 
reduced for solution of the spin axis by the CFLS method. No bias corrections were 
applied, so that comparison is made with the uncorrected, zero bias IDCT results. Also, 

14 orbits of SAS-2 were selected for which star sensor and IDCT solutions exist, but the 
best evaluations of bias corrections were applied before making the CFLS solutions. 

Star Sensor 

The Attitude Determination and Control Section produced star sensor solutions in a real- 
time, operational environment necessitating fast evaluations of the data. While the precision 
*Meyers, G. F. et aL, SAS-1 Attitude Data Analyses, GSFC X-542-71-363, 1971. 
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of the instruments to one arc-minute indicates accurate solutions, a comparison with 
definitive star sensor solutions for a few orbits was made to detect any unsuspected varia- 
tions in the accuracy. This comparison showed the operational solutions to be within a 
few hundredths of a degree of the definitive star sensor solutions (table 5). This accuracy 
is sufficient for operational star-sensor spin-axis solutions, to serve as a base for determining 
the accuracy of a cone intersection solution by the IDCT or CFLS methods. 


Table 5 

Comparison of Operational with Definitive Star Sensor Solutions 


Orbit 

; Number 

1 

Operational Spin Axis 

Average* Definitive Spin Axis 

Arc 

Difference 

Right 

Ascension 

Declination 

Right 

Ascension 

Declination 

386 

270?423 

-25f 256 

270? 340 

-25?314 

Of 100 

416 

273?697 

-16?336 

273?674 

-1 6f 458 

Of 123 

476 

279?792 

- 1 f 280 

219° 61 A 

- 1 f 333 

0°052 

491 

28 1 f 829 

- 1? 102 

281?836 

- If 100 

0?003 

523 

3 1 2f 3 1 0 

5f 891 

3 1 2f 359 

5°709 

0? 189 

537 

306?664 

6?615 

306?643 

6? 5 80 

0?034 

584 

336?390 

-1 2f 590 

336?307 

-1 2f 67 1 

Of 1 1 5 

599 

339?533 

-20?379 

339?472 

-20?470 

0? 108 


* Average values of right ascension and declination are given. 
Deviation in right ascension was between ±0?059 and ±0?182. 
Deviation in declination was between ±0?048 and ±0?158. 


Magnetometer-Sun Sensor 

To satisfy mission requirements the magnetometer-sun sensor cone angle solutions of the 
spin axis must be within 5° of the actual, definitive results. To obtain this accuracy by the 
IDCT, it is necessary to use magnetometer and sun data selected throughout an orbit 
between changes of the spin-axis orientation. A similar but slightly different procedure is 
followed for CFLS solutions. Table 6 gives the IDCT and star solutions by orbit for the 
28 playback orbits which were analyzed. There are three solutions under IDCT representing 
the three different residual bias considerations. Specifically, these are (1) the case with no 
residual magnetic field correction; (2) the case with no bias in the x and y magnetometers 
and a 239 milliampere/meter (+3 millioersted) bias in the z-magnetometer; and (3) the case 
with no bias in the x, a +3 millioersted (moe) bias in the y, and a 398 milliampere/meter 
(+5 millioersted) bias in the z-magnetometer. The no-bias corrections case, (1) above, 
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Table 6 

Comparison of Operational Star Solutions with IDCT Solutions for Different 

Residual Magnetic Field Biases 
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produces adequate solutions; the (0, 0, 3) millioersted corrections case, (2) above, corres- 
ponds to the prelaunch determination and improves the solution; but the in-flight deter- 
mined (0, 3, 5) case, (3) above, clearly is the best with a mean arc error of 0?47 compared 
to 1?50 and 0?81, respectively. 

Table 7 is the set of CFLS solutions compared with star solutions for 1 1 of the 28 orbits 
for which raw telemetry and orbit data were specifically reduced. The distinction between 
the two different CFLS solution columns is that for the CFLS/PQ case the input equations 
consisted of the sensor cone angles and the vector components of the sun and magnetic field 
only, but for the PQV case the normal vector V and its associated cone angle were added to 
the equations. From this table the CFLS/PQV approach is seen to be much stronger. The 
reason for this is that the CFLS/PQ equations do not produce good geometrical constraints 
on the solution, since the magnetic field vectors are all within a few degrees of arc and the 
sun vectors vary by only one degree. This condition magnifies errors and creates the inac- 
curate CFLS/PQ solutions. The orbits 386, 387, 553, and 554 had data selected from 
throughout the orbit and, therefore, had a wider range of magnetic field vector positions, 
which resulted in significantly more accurate solutions. The CFLS/PQV equations provided 
strong geometrical constraints and resulted in significantly more accurate solutions. 

Table 7 

Comparison of Operational Star Solutions with CFLS Solutions 



Star 

CFLS/PQ 

CFLS/PQV 

| Orbit 

Right 

Ascension 

Declination 

Right 

Ascension 

Declination 

Right 

Ascension 

Declination 

370 

327?78 

-30? 01 

33 8? 54 

-1 8f 57 

328?68 

-28?76 

386 

270. 75 

-25. 25 

272. 53 

-23.09 

270. 83 

-23.25 

387 

270. 83 

-25.25 

262. 96 

-23.42 

271.30 

-22. 42 

416 

274. 00 

-16.32 

280. 85 

-12.33 

275. 12 

-13.73 

432 

268.05 

-12.78 

273. 94 

-10.38 

269. 84 

-10.78 

446 

275. 64 

-11. 13 

248. 00 

-14.58 

279. 36 

- 8.02 

447 

275.57 

-11.09 

272. 36 

-10.39 

275. 80 

- 9.85 

553 

300. 23 

8.81 

299. 02 

10. 20 

296.33 

9. 14 

554 

300. 20 

8.74 

297. 28 

9.81 

295. 24 

9. 05 

577 

337.02 

-12. 83 

348.31 

-14. 46 

336.52 

-11.71 

599 

339. 82 

-20. 26 

340. 29 

-13.26 

341.43 

-19. 18 


Comparison of SAS-1 IDCT and CFLS Solutions 

For direct comparison of the IDCT with the CFLS results, table 8 is given. Qualifying 
explanations of the CFLS solutions are necessary, however, to interpret the results. From 
tables 6,7, and 8, it is apparent that the CFLS solutions are significantly less accurate than 
even the no-bias IDCT values, but operational conditions account for some of the variation. 
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Table 8 

Comparison of CFLS with IDCT Solutions 



Star 

CFLS/PQV 

IDCT (0,0,0) Bias 

Orbit 

Right 


Right 


Right 



Ascension 

Declination 

Ascension 

Declination 

Ascension 

Declination 

370 

327? 78 

-30? 0 1 

328?68 

-28?76 

328?39 

-29?06 

386 

270. 75 

-25.25 

270. 83 

-23. 25 

270.51 

-24. 32 

387 

270. 83 

-25.25 

271.30 

-22.42 

270. 87 

-24. 42 

416 

274. 00 

-16.32 

275. 12 

-13.73 

274. 67 

-15.20 

432 

268. 06 

-12.78 

269. 84 

-10.78 

268. 58 

-11.53 

446 

275.64 

-11. 13 

279. 36 

- 8.02 

276. 19 

-10.22 

447 

275.57 

-11.09 

275. 80 

- 9.85 

275.54 

-10. 09 

553 

300. 23 

8.81 

296. 33 

9. 14 

299. 59 

9.38 

554 

300. 20 

8. 74 

295. 24 

9.05 

299. 89 

9. 50 

577 

337.02 

-12.83 

336. 52 

-11.71 

336. 96 

-11.93 

599 

339.82 

-20. 26 

341.43 

-19. 18 

339. 75 

-19.21 


That is- 

• For all orbits except 386, 387, 553, and 554, data were selected from a few short 
time intervals rather than scattered evenly throughout the full orbits. The short time 
intervals were about 5 minutes in duration and corresponded to periods of properly 
correlated time between universal time and the spacecraft clock. 

« The sun angle readings were accurate only to within l°;yet, alterations of Of 1 in 
tests of orbits 553 and 554 produced deviations of arc from 1° to 2°. Daily records 
and analysis would allow interpolation of the sun angle to an accuracy of less than 
0? 5 . (For IMP-6 such interpolation provides sun angle values good to 0? 05.) 

• The magnetic field magnitudes computed from telemetry data differed from the field 
magnitudes supplied by the Orbit Determination and Control Section by up to +1600 
to -800 milliamperes/meter (+20 to -10 millioersteds). The field quantities derived 
from telemetry were used in computing cone angles for consistency and to minimize 
any systematic error; but as has been shown in the residual bias comparison on the 
IDCT solutions, a remnant difference of a few millioersteds affects the solutions 
proportionately. 

• Interpolation from telemetry time to the nearest minute was not made on the meas- 
ured magnetic field components to bring them to the same time as the magnetic field 
model inertial coordinates, because the measured z-components were not varying 
significantly in those few seconds; but the x- and y-components, which had to be 
interpolated, varied to a greater degree than the field model would indicate as likely. 
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Table 9 


CFLS Solution Errors in Right Ascension, Declination, and Arc 
and the Standard Deviation in Arc for Each Orbit 


Orbit 

ARight 

Ascension 

ADeclination 

AArc 

(T Arc 

370 

+0?89 

+ If 25 

1 f 47 

0?72 

386 

+0.08 

+2.00 

2.01 

0. 74 

387 

+0. 47 

+2.73 

2.86 

0,71 

416 

+1. 12 

+2. 59 

2.81 

0. 83 

432 

+1.78 

+2.00 

2.66 

0. 72 

446 

+3.72 

+3. 12 

4.82 

0. 94 

447 

+0. 23 

+ 1.24 

1.26 

0. 94 

553 

-3.90 

-0. 23 

3.87 

0. 745 

554 

-4. 96 

-0.31 

4.91 

0. 70 

577 

-0.50 

+ 1. 12 

1.22 

0. 945 

599 

+ 1.61 

+ 1.08 

1.86 

0. 90 


Error Analysis 

Table 9 is a summary of all of the errors in the CFLS/PQV solutions. Typical errors are 
given by the differences between the CFLS solutions and the star solutions in right ascen- 
sion, declination, and arc. 'These differences are equivalent to the net effect of the errors 
caused by individual data point contributions. The means of these errors are as follows: 


A RA 

= -HD?02 

A Dec 

= +1 f 5 1 

A Arc 

= 2°1Q 


Note that the trend of these errors is positive, indicating a systematic error presence con- 
sistent with the comparison qualifications described in the section above. 

The strongly negative differences for orbits 553 and 554 indicate the possibility of another 
factor at work in these orbits, and this in fact is the case. Since data for 553 and 554 were 
selected at several different time intervals throughout the orbit, the solutions for these orbits 
were expected to be superior. Checking the CFLS/PQ solutions (table 7) shows that the 
best results for 553 and 554 were obtained by that procedure: 


Orbit: 

553 

554 

Right ascension: 

- 1 f 2 1 

- 2?92 

Declination: 

+1?39 

+0?77 

Arc: 

1?84 

3°00 
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Examination of the CFLS/PQV equations shows that the geometry for the computation of 
the V-vector-to-spin-axis cone angle involved the sum of two spherical triangle legs ((3 and 
S), which was less than the third leg ( 7 ?). This condition existed due to subtractive errors 
in the measurements of (3 and 5, combined with the situation of the spin axis becoming 
aligned with the plane of the sun and the magnetic field vector, making |3 + 5 « rj, with 
large variations in the cosines of small angles for small deviations in the angles. An alter- 
native angle determination scheme can be incorporated in the algorithm to bypass this 
problem. 

The root-mean-square (rms) error can be determined for any data which has a mean, and it 
is useful for finding the range about an arbitrary origin in which new data is likely to fall. 
For the situation where the data is taken as the spin-axis difference components between 
the star solution and the cone angle solution for each orbit, the rms is 

(4 ‘ 2)/N 

where A. = one of the difference components ARA., ADeCj, or AArc p and 
N = the number of total orbits.* 

These values are 


™s ARA 

= 2°33 

RM S ADec 

= 1°85 

^AArc 

= 2? 99 


and the ARA have been normalized to degrees of arc. (See the section on SAS-2 error 
analysis for further discussion.) 


The standard deviations of the angle components for the entire set of orbit solutions takes 
the form— 



(A i - a > 2 


N - 1 


where A. is as defined above, A is the average of that component over all orbits, and N-l is 
is the degree of freedom (one constraint being used in the determination of the mean A). 


*Kirby, P. et al., Data Processing Plan for Small Astronomy Satellite (SAS-Aj, GSFC X-564-71-130, December 1970. 
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The properties of the standard deviation are well known, such that the selection of a new 
orbit data set should produce a solution which differs from the appropriate star solution in 
right ascension by less than a ARA , in declination by less than a ADec , and in arc by less than 
a AArc> 68 percent of the time in each case. These values are— 


a ARA 


2°43 


a ADec 


If 11 


An indirect computation of the standard deviations of the parameters from a least squares 
solution can be obtained without resorting to separate spin-axis solutions of each set of 
sensor measurements. These standard deviations are derived from the inverse of the normal 
matrix, N' 1 (= [U T U] * 1 ; see the section above, “The Closed-form Least Squares Method” 
and Appendix D). This matrix is called the variance-covariance matrix of the solution. The 
diagonal elements of N' 1 are extracted, and are the normalized, nondimensional variances 
of the solution parameters. To obtain the required variances (squares of the standard devia- 
tions), a multiplicative factor— the square of the standard error of a single measurement- 
must be introduced. If d^ , d^ , and d^ are the diagonal elements, then the true variances 
are 


'a 2 = m 2 d 2 ; 'a 2 = m 2 d 2 ; o' 2 = m 2 d 2 

xx xx 9 yy yy * zz zz 


If the variances are desired in spherical rather than Cartesian coordinates, differentiation 
properties apply such that 


RA 


± ( Xg yy - y g xx) 
x 2 + y 2 


and 


Dec 


(1 - z 2 ) J 


2\V2 


where 


and 


where 


f 


from o (tan a) = a (x/y) = 


cos 2 a 



cos 2 a = 


x 2 + y 2 


o (sin 6) = a (z) = cos 5 a 


cos 


5 = \/\ - Z 2 


]■ 
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However, a more useful variance is the net or total arc variance, 

°a ~ (p 2 + o 2 + a 2 ), 

Arc v xx yy zz 

which, with the appropriate multiplicative factor, gives the distribution of angular error 
around the spin axis solution. That is, the standard deviation 


is the radius around the solution of the spin-axis vector for a particular orbit, within which 
another spin-axis solution, based on similar measurements, would fall 68 percent of the 
time. A good estimate of m for the determination of a Arc is the average of the magnetome- 
ter standard error and the sun sensor standard error— (If 5 + 0?5)/2 = 1°— for all sensor 
readings in each orbit. The results of this assumption are shown in table 9 for o Arc . 

Peak arc errors of about 5° for orbit 554 and about 4° for 553 were examined above and 
were attributed to a geometrical configuration of the sun, magnetic field, and spin axis vec- 
tors, which, under the present algorithm, tend to propagate errors in this situation rather 
than damp them out. 

The peak arc error of nearly 5° for orbit 446 is clearly related to the CFLS/PQ solution 
arc error for that orbit of nearly 30°. Obviously, a gross error is present in the data, though 
a careful check of all values indicates that whatever error is present is due to inaccuracies 
in the raw data. 

Continued analysis and refinement of the CFLS method as applied to SAS, and increased 
care in the solution and reduction of data, would reduce error peaks to the present mean 
value, and would tend to produce solution accuracies approaching the IDCT results. Five- 
degree accuracy is obviously obtained, and operational accuracies of 2° or better are indi- 
cated and were achieved for SAS-2 data, as will be explained below. 

SAS-2 Attitude Solutions 

The SAS-1 CFLS spin-axis solutions in the section, “SAS-1 Attitude Solutions,” were based 
on certain prelaunch determined spacecraft parameters, and were discussed with no consid- 
eration of bias or calibration adjustments to the magnetometer data. This approach is 
reasonably accurate, but incorporation of bias and recalibration values significantly improve 
the solutions. The results obtained by using such computed adjustments will be considered 
in this section. 

SAS-2 spin axis solutions for certain orbits for which IDCT solutions exist were redetermined 
using the CFLS method. In addition to the no-bias and nominal calibration solutions, CFLS 
solutions based on residual spacecraft magnetism and sensor biases and new millivolts-to- 
millioersted calibration equations (as determined by the Attitude Determination and Control 
Section) were obtained. Comparison is made here, for each orbit, of the operational star 
solutions (considered as accurate) with the nominal calibration (discussed in the section, 
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“IDCT Solution Applied to SAS-1 ”), the no-bias correction IDCT and CFLS solutions, and 
with the best bias, recalibrated CFLS solutions. 

IDCT solutions with the best bias and recalibration values were not available for the selected 
orbit data, but another group of orbits for which such IDCT solutions exist provided statis- 
tics for comparison with CFLS solution statistics. Nonorthogonality and misalignment 
adjustments were found to be very small and were ignored. 

Bias and Recalibration Parameters 

Determination of SAS-2 residual magnetic field bias corrections and millivolts-to-millioersted 
(milliampere/meter) calibration parameters is made by solving a single equation for each axis. 
The equation is 

f.V. = sin 0. cos (cot + <t > 0 ) (H x sin 5 cos a + H y sin 5 sin a - H z cos 6) 

+ sin 0j sin (cot + 0 o )(-H x sina + H y cos a) 

+ cos 0j (H x cos 5 cos a + H y cos 8 sin a + H z sin 5) + B. 

where H x , H y , and H z = the inertial components of the earth’s magnetic field, deter- 
mined at the spacecraft position from a theoretical model of 
the earth’s magnetic field; 

V. = the voltage in millivolts from a particular magnetometer axis 

(i = x,y,z); 

a and 5 = the right ascension and declination components, respectively, 

of the spin axis; 

co = the spin rate; 

0. = the alignment angle of the ith magnetometer axis; 

</> o = the initial (t = 0) azimuth of the ith magnetometer axis; 

i 

f. = the slope of the ith magnetometer axis calibration curve (in 

moe/mv), and 

B. = the ith magnetometer axis bias (in moe). 

The form of this equation results from transforming the inertial coordinates of the magnetic 
field model vector into the spacecraft reference frame (spin axis = z axis), and projecting 
these components onto each of the magnetometer axes. (See figures 18 and 19.) This is 
equivalent to taking the dot product between the magnetic field vector in spacecraft coor- 
dinates and the magnetometer axis unit vector. The B term is the linear addition of the 
residual bias. 
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Figure 18. Magnetic Field Vector H and Spin Figure 19. Magnetometer Axis M in Spacecraft 

Axis in S in Inertial Coordinates Coordinates 

The equation can be solved in a least squares sense for any of the last seven parameters, 
since over a given time interval within an orbit they are nearly constant and the equation 
can be overdetermined. 

For final determination of f and B. (i = x, y, z), selection of the best values for the other 
parameters from prior solutions and other methods is made, and an iterative, differential 
correction least squares fit is performed on the equations with the other parameters held 
constant. Results were taken over many orbits and averaged to the nearest millioersted for 
B. and to the moe/mv ratio of whole numbers nearest the nominal ratio of 700/500 for f.. 

The initial determinations were for the biases only— not including the recalibration slopes— 
and produced the residual magnetic field biases of x = 0, y = 0, and z = -7. While the mag- 
nitudes of the magnetic field computed using these biases were better than fields computed 
using the (0, 0, 0) bias set, when compared with the field model values, the spin axis solu- 
tions derived from the (0, 0, -7) set were not as close to operational star solutions as were 
solutions based on the (0, 0, 0) set. This was interpreted as a requirement for computation 
of the calibration slopes. The resultant recalibrated millivolts-to-millioersted conversion 
equation for the x- and y-axes was 

millioersted = (millivolts + 254) (728/508) - 364 


and for the z-axis was 

millioersted = (millivolts + 254) (764/508) - 382. 

The residual bias corrections were x = -1 , y = +1 , and z = +2. This set of bias corrections 
and calibration slopes led to more accurate attitude solutions relative to the star solutions 
and to magnetic field magnitudes closest to the field model. 
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The presence of bias adjustment values, however, does not necessarily indicate an actual 
residual magnetism on the spacecraft, but could mean that the millivolts-to-millioersted 
conversion calibration is offset such that a zero millivolt reading actually reflects a nonzero 
field value. 

Comparison of SAS-2 Attitude Solutions by Various Methods 

Operational star, IDCT, and CFLS solutions are given in table 10 for 13 orbits. The IDCT 
solutions are based on the nominal calibration (slope of curve = 700/500) and no residual 
bias computation of the magnetic field from the magnetometer readings. The CFLS solu- 
tions are based on both the nominal calibration (no-bias computation and a nominal cali- 
bration, z-axis bias of -7) and a recalibrated (biased by -1 , 1 , 2) field computation, and are 
listed separately. The CFLS solutions are all arrived at by the CFLS/PQV approach dis- 
cussed in the section on the magnetometer-sun sensor. 

The arc differences show a clear improvement in solutions derived from the recalibrated, 
biased CFLS approach over the nominal calibration cases. Solutions with the (0, 0, -7) bias 
are less accurate for nearly every orbit, yet the magnetic field vector magnitudes computed 
with this bias are closer to the field model vector magnitudes than the magnitudes of the 
(0, 0, 0) biased vectors. The average magnitude differences over each orbit between mag- 
netometer-measured field vectors and the corresponding field-model vectors are given in 
table 1 1 . 

Table 1 1 reflects the area of improvement of solutions when a jointly derived bias set and 
recalibration are included in determination of magnetometer-measured fields. A slight reduc- 
tion of the -1 , +1 , +2 biases causes the “magnetometer-model” differences to be more posi- 
tive, hence closer to zero, and the magnetometer-derived magnetic field magnitudes would 
be closer to the model values. Improvement of the spin axis solutions, however, is not 
assured— compare the (0, 0, 0) and (0, 0, -7) results— and field magnitude differences over 
the entire range of 2000 orbits for which the bias/recalibration determination was made 
should show both positive and negative values, and the average should tend to zero. Other 
parameters possibly contributing to attitude solution errors are considered more fully in 
the next section. 

For comparison with the recalibrated, bias-derived CFLS solutions, a separate set of orbits 
solved using the recalibrated, bias-derived IDCT is necessary, since such solutions for the 


45 



Table 10 

Comparison of Operational Star Solutions with IDCT and CFLS Solutions 





















Table 1 1 

Comparison of Magnitudes of Magnetic Field Model Vectors 
with Magnetometer-derived Field Vectors 


Orbit 

Averages of Differences Between Magnetometer-derived 
and Field Model Magnitudes 

(0, 0, 0) Bias 
Nominal Calibration 

(0, 0, -7) Bias 
Nominal Calibration 

(-1, +1, +2) Bias 
Recalibration 


moe* 

mA/m 

moe 

mA/m 

moe 

mA/m 

1063 

- 7.68 

-611.15 

-5.67 

-451.20 

-0.37 

- 15.31 

1068 

-7.77 

-618.32 

-4.98 

-396.29 

-0.52 

- 41.38 

1072 

-9.39 

-747.23 

-6.01 

-315.13 

-0.59 

- 46.95 

1083 

-8.77 

-697.89 

-6.33 

-503.72 

-1.08 

- 85.94 

1088 

- 8.75 

-696.30 

-5.56 

-442.45 

-0.04 

- >3.18 

1093 

-8.89 

-707.44 

-6.31 

-502.13 

-0.84 

- 66.84 

1198 

-9.15 

-728.13 

-6.25 

-497.36 

-1.54 

-122.55 

1103/1104 

-9.27 

-737.68 

-6.00 

-477.46 

-0.83 

- 66.05 

1108 

-9.21 

-732.91 

-5.65 

-449.61 

-0.57 

- 45.36 

1113 

-8.01 

-637.41 

-5.35 

-425.74 

-0.55 

- 43.77 

1118 

-9.27 

-737.68 

-6.00 

-477.46 

-0.69 

- 54.91 

1124 

-9.14 

-727.34 

-6.57 

-522.82 

-1.39 

-110.61 

1128 

-7.80 

-620.70 

-5.10 

-405.84 

-0.75 

- 59.68 


♦Millioersteds 


orbits of table 10 were not determined. Table 12 reproduces these results, and an analysis 
of them with respect to the data of table 10 is provided in the following section. 

Error Analysis 

Deviations of the magnetometer/sun sensor spin-axis solutions from the operational star- 
sensor solutions are expected because of uncertainties and errors in the sensor measurements. 
Random errors tend to cancel when many measurements are used, and systematic errors 
lessen as the uncertainties in the measurements are resolved. Improvement of solutions due 
to recalibration and biasing of the magnetometers has been noted above. 
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Table 12 

Recalibration and Bias-derived IDCT Solutions 
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Further improvement of solutions is expected if sun angle values are refined. This refine- 
ment is possible by reducing the systematic error in the sun angle measurements created by 
the ±0?5 resolution limitation of the sun sensor. That is, a time history of the sun angle 
measurements and interpolation between the sun angle changes can produce values good to 
0 ? 1 . 

With the sensor measurements corrected or specified to the precision limits of the sensor, 
a residue systematic error and any net random error remain. Table 13 presents an error 
summary of solutions by the various methods given in the previous section, compared with 
operational star solutions. All of the measures of dispersion indicated in table 1 3 relate to 
the components of the solution differences in such a way as to point up systematic errors, 
inaccuracies, and imprecisions in the measurements involved in the solutions. 

The averages are the arithmetic means of the differences between the sun sensor-magnetom- 
eter and star sensor solution components (right ascension, declination, and arc separation), 
and indicate the accuracy of the set of cone angle solutions with respect to the set of star 
solutions. Large (>0f 5) dispersions in the right ascension and declination components also 
reflect a probable systematic error. 

The rms deviations are measures of dispersion of the weighted magnitudes of the differences 
and indicate the range about the star solution component within which an arbitrary cone- 
angle-solution component is most probable. This is another estimator of the accuracy of 
the cone-angle solutions. 

The standard deviations are measures of dispersion of the solution differences about the 
mean of the differences, and thus serve as an indicator of the precision of the solutions. 

If the average c omponent difference is zero, then the a and rms entries coincide except 
for a>/N/(N-l) term which must be factored from o (where N is the number of solutions, 
N= 13 .-Vn/(N-1) =^13712= 1.04). 

The arc separation and declination differences are consistent statistics for any solutions, in 
that a degree of arc or declination is the same regardless of coordinates. A degree of right 
ascension, however, varies in size in proportion to the declination. To cancel this effect, 
right ascension error statistics have been normalized and are expressed in degrees of arc for 
valid comparison. 

The ARA averages are close to zero for all methods, with the ARA’s rms and a entries 
accordingly almost coinciding. However, the rms and a magnitudes from the CFLS methods 
are around 0?4, signifying a large range in ARA values. 

The ADec averages show considerable variation, with the recalibrated, biased CFLS solution 
average, being significantly smaller than the others. In addition, the rms and a entries for 
that method are small, and show high accuracy and precision due to proper data adjustment. 

The large differences between the rms and a values for the other CFLS methods, with a 
being much better than rms, reflect the systematic errors present in the improperly or 
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Table 13 

Statistical Comparison of the Accuracies of IDCT and CFLS Solutions with Respect to 

Operational Star Solutions 



50 


ARA = Averages of the differences in right ascension between the star and the magnetometer-sun sensor solutions. 

RMS Dev. of ARA = Root-mean-square deviations in the right ascension differences. 

a ARA = Standard deviations of the right ascension differences about the mean of the differences. 

Based on the separate set of orbits of table 12. 























unadjusted approaches. The IDCT declination statistics are good, but seem significantly 
less accurate and precise. 

The arc statistics, on the other hand, suggest that the IDCT method gives slightly better 
solutions than the best CFLS approach, but the average arc deviation and a ±1 standard- 
deviation range around the average for each method shows the IDCT range to be encom- 
passed by the range from the biased, recalibrated CFLS solutions. Any preference of 
method, then, seems to be unwarranted, based on these statistics from nonidentical orbits. 

A few additional statements should be made to further interpret or clarify the statistics. 

• An error in the sun angle has its strongest effect on the right ascension of these spin 
axis solutions, due to the constraint that the sun’s declination is always within ±23?5, 
and to the fact that the spin axis declination is approximately -21° . Thus, an inaccu- 
racy in the sun angle expands or contracts the circle of possible solution points, which 
in the neighborhood of these solutions is primarily a change in right ascension. 

Similarly, the direction of a magnetic field vector detected from the SAS-2 orbit is 
constrained between +60° and +80° declination, which causes any magnetic field 
error to show up primarily as an error in the calculated spin-axis declination. 

From these conditions, it is apparent why the solution statistics improved so signifi- 
cantly in the declination but not in the right ascension of the biased, recalibrated 
CFLS solutions. The significantly better values in the right ascension statistics for 
the IDCT method possibly can be attributed to two factors: (1) Different sun posi- 
tion-spin axis relationships lessen the right ascension sensitivity to sun angle error in 
the IDCT data. (2) Sun angle interpolation may have been employed on the IDCT 
data to give better than ±0?5 accuracy. This was not done with the CFLS reduced 
data. 

• The arc is a magnitude measurement with only positive contributions to arc statistics; 
thus no cancellation of positive and negative elements is possible. This means that 
random effects cannot cancel out, but must be a part of the final net averages. Simi- 
larly, the right ascension and declination rms and a statistics have no cancellations 
(since terms are squared), and so include all random error components. 

• The magnetometer and sun sensor data were selected independently of the star sensor 
data for each orbit; therefore solutions based on perfect magnetometer-sun sensor 
readings and calibrations would not be identical to star sensor solutions. This is 
because of the slow drift of the spin axis vector over an orbital period, and the resul- 
tant, slightly-differing solutions based on data taken at different times. 

• The standard deviations of the spin-axis spherical coordinates of each orbit were 
computed directly from a series of biased, recalibrated PQV solutions for each set of 
sensor measurements related to the CFLS (-1, 1, 2) solution, in contrast to the 
variance-covariance matrix approach of the section on SAS-1 error analysis. These 
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I statistics are presented in table 14. Note the small variations in o RA . This is a result of the 
strong dependence of the RA solution components on the sun sensor measurements, 
and since the sun sensor showed no change over a single 95-minute orbit, a high pre- 
cision in this component is expected. 

The o" Arc and o Dec elements accordingly almost coincide. These statistics are reasonably 
precise and give evidence that the magnetometer measurements as a group were at the 
limit of resolution, with few, if any, bad readings. 


Table 14 

Orbit Statistics from Biased, Recalibrated PQV Statistics 


Orbit 

°Aic 

°RA 

°Dec 

1063 

0?39 

0?02 

0?39 - 

1068 

0.50 

0.04 

0.50 

1072 

0.50 

0. 04 

0. 50 

1083 

0. 59 

0. 03 

0.59 

1088 

0. 48 

0.04 

0. 48 

1093 

0. 29 

0. 03 

0. 29 

1098 

0.59 

0. 06 

0.59 

1103/1104 

0. 80 

0.07 

0. 80 

1108 

0.45 

0.04 

0. 45 

1113 

0.68 

0. 06 

0. 68 

1118 

0. 48 

0. 03 

0. 48 

1124 

0.56 

0. 05 

0.56 

1128 

0. 77 

0.05 

0.77 


Advantages and Disadvantages of the Four Approaches 

The four approaches to cone intersection solutions of the spacecraft spin axis presented in 
this report are all in use on one or more of the satellite data reduction systems at Goddard 
Space Flight Center, and each has limitations or liabilities as well as the proven capability 
of producing satisfactory solutions. 

The critical problem for the cone intersection methods for spacecraft without the time 
measured between sensor readings occurs as the two intersection points of two cones 
approach each other and the initial estimate is not necessarily closest to the true solution. 
For the IDCT and CFLS/PQ least-squares approaches, data separated in time as widely as 
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feasible provide multiple vectors, with sufficiently different cone angles, to select unam- 
biguously the proper solution. The non-least squares approaches, producing multiple 
solutions for each complete telemetry time! interval based on each ambiguity in the data, i 
must be augmented by time histories showing interval-to-interval status of the solutions. 

In this way, however, slow changes in the spin axis due to external torques can be detected 
and followed. 

In order for the least squares approaches to follow a slowly-moving spin axis, alterations in 
the algorithms have been designed which utilize the overdetermined nature of the equations 
to solve for the initial position, RA„ , Dec n , or X n , Y n , Z n ; the rates of change, RX, Dec, 
or X, Y, Z; and rates of change of the rates of change, RA, Dec, or X, Y, Z; such that, for 
some time At after the initial time t 0 , the spin axis position components are given by 

RA = RA 0 + R'A • At + RA • At 2 
Dec = Dec Q + Dec • At + Dec • At 2 
or 

X = X 0 + X • At + X - At 2 

Y = Y 0 + Y • At + Y • At 2 

Z = Z Q + Z • At + Z • At 2 . 

The greatest limitation on the least-squares approaches occurs during periods of rapid spin- 
axis reorientation maneuvers, when the overdetermined telemetry data is inconsistent with 
good attitude solutions. 

The chief difference in the least-squares methods is in the computation complexity of the 
IDCT over the CFLS. By relying on a single pass through the algorithm steps, the CFLS 
method is faster than the more complex and iterated IDCT, although solutions may be 
slightly more accurate from the IDCT. 


Goddard Space Flight Center 

National Aeronautics and Space Administration 
Greenbelt, Maryland July 1974 
311-07-14-01-51 
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APPENDIX A 


EULER ANGLES 


A spacecraft with body-fixed axes xyz and an inertial system with coordinates XYZ is con- 
sidered. (See figure A-l .) To transform xyz into XYZ, a series of three rotations about 

initial and intermediate axes is the maximum 
number required, regardless of the relative 
orientation of the two sets of axes. These 
rotations are known as Euler angles. They 
are not unique. That is, the general solution 
of the transformation of any right-handed 
orthogonal coordinate system into another 
coordinate system can be characterized by 
several different sets of Euler angles. The 
first rotation can be about any axis— x, y, or 
z. If the first rotation axis selected is the z, 
then several angles of rotation leading to the 
desired transformation are possible. Specifi- 
cally, these angles may be the rotations of y 
into either the ZY, zY, zX, ZX, or zZ planes. 
There are two intersections of the rotation 
of y about z with each of the above planes 
giving ten possible choices of the rotation. Similarly, there are 10 possible choices for 
the rotation of the x-axis about the z-axis into the five planes. 

This first rotation is necessary in order to position the transformed axes with respect to 
the fixed axes, in such a way that the second rotation brings one transformed axis into 
coincidence with the desired fixed axis. The first rotations indicated above satisfy this 
condition for one of two reasons. Either the rotation of an axis into a plane formed by 
a fixed axis and the axis of rotation places the third transformed axis perpendicular to 
that plane from which the second rotation about the perpendicular axis is made; or the 
rotation of an axis into a plane formed by two fixed axes allows the second rotation to 
be about this axis, such that the appropriate axis is made both perpendicular to the fixed 
axis plane and coincident with the third fixed axis. The second rotation having trans- 
formed one axis into one fixed axis, the third rotation is about these coincident axes 
and completes the transformation. 

Obviously, if the first axis of rotation had been the x or y rather than the z, parallel 
statements could be made. Thus, 30 first rotations are possible. The second and third 
rotations are dictated by the first, resulting in the 30 different transformations. 


z 



Figure A-1. Spacecraft xyz Axes and Inertial 
Coordinates XYZ 
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As illustrations of the above (see figure A-2), the following rotations are considered: The 
transformation of x-, y-, z-axes into X-, Y-, Z-axes by an initial rotation around the z-axis 
(by angle 0), moving the x-axis into the XY plane; second, a second rotation about the 
new x-axis until the z-axis coincides with the Z-axis (angle 0); and third, a final rotation 
around the z=Z axis (by an angle \p). 



(a) Initial Rotation about 
z-axis 


(b) Rotation about new 
x'-axis 


(c) Final Rotation about 
Z-axis 


Figure A-2. Rotation of Spacecraft Axes 


Each rotation can be described mathematically by analytic geometry or by matrix theory. 
Thus, 


z' = z; z" = (-sin0)y' + (cos0)z'; 

x = (cos0)x + (sin0)y; x"= x'; 

y' = (-sin 0) x + (cos 0) y; y"= (cos 0) y' + (sin 0) z'; 


x'" = (cos 0) x" + (sin 0) y" = X 
y'" = (-sin 0) x" + (cos 0) y" = Y 


or more compactly, 
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The total transformation is even more convenient by using matrix notation. Instead of 
algebraically substituting to get 

Z = Z (x, y, z, 0, 0, 0); X = X (x, y, z, 0, 0, \]j); Y = Y (x, y, z, 0, 0, i p), 


the matrix product- 



cos 0 sin 0 0 

-sin 0 cos 0 0 

0 0 1 


10 0 cos 0 sin 0 0 

0 cos 0 sin 0 -sin 0 cos 0 0 

0 -sin 0 cos 0 0 0 1 




cos 0 cos 0 -sin 0 cos 0 sin 0 cos 0 sin 0 + sin 0 cos 0 cos 0 sin 0 sin 0 

-sin 0 cos 0 - cos 0 cos 0 sin 0 -sin 0 sin 0 + cos 0 cos 0 cos 0 cos 0 sin 0 

sin 0 sin 0 -sin 0 cos 0 cos 0 



provides the desired relationships. These are the Euler rotations used in Goldstein 
(reference 1 , p. 107), with constraints 0 < 0 < 27r, 0 < 0 < 7r, 0 < 6 < 2n. 

An important convention should be noted. Angles are considered positive for counter- 
clockwise rotations and negative for clockwise rotations. The rotations in illustration A-2 
are ccw for 6, and cw for 0 and 0; therefore, the angles 0 and 0 are negative. 

Next, the following rotations are considered (figure A-3): The transformation of xyz into 
XYZ by an initial rotation around the x-axis (by an angle £), moving the z-axis into the xZ 
plane; a second rotation about the new y-axis (perpendicular to the xz'Z plane), until the 
z'-axis coincides with the Z (angle f); and a final rotation around the z" = Z axis (by an 
angle i?), which brings x into X and y into Y. 
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APPENDIX B 
SPACECRAFT SENSORS 


SUN SENSORS 

The sensor angle derived from the digital sun sensor is the angle between the spacecraft spin 
axis and the vector from the spacecraft to the sun. This vector is computed from ground 
information of the instantaneous positions of the spacecraft and the sun in geocentric iner- 
tial coordinates. 

The most commonly used digital sun sensor has photocells arranged in 10 columns, pro- 
duces nine bits, and indicates the sun angle to within 0?5. To minimize the number of 
photocells, the array design places the cells in a pattern which requires various-sized cells 
overlapping many rows. One such pattern is shown in figure B-l . 

The long dimension of the box-shaped sensor is mounted perpendicular to the spin axis 
and is centered in the xy-plane. There are two slits with 180° fields of view, such that 
every point on the celestial sphere is within the fields of view once each spin period. A 
long sun-angle slit, in the xy-plane, permits the sun’s rays to fall on a particular row of cells 
across nine columns, and the array is constructed in such a way that the resulting 1- and 
0-bits indicate a particular sun angle. If a photocell in a particular column is excited, the 
corresponding binary digit in a telemetry word is a 1-bit. If the light passes over a blank 
cell, the corresponding bit is a 0-bit. The bits are read out in order by the telemetry elec- 
tronics, and eventually are decoded on the ground. 

The second slit, the command slit, is parallel to the spin axis. The sun image from this slit 
falls only on the first column of photocells. When the sun image first falls on this column 
an electronic pulse initiates a counter, and the next sequential pulse stops the counter. The 
sensor electronics centers each pulse with respect to the total time that the sun is viewed by 
the particular photocells, and the resulting total count between two pulses provides data on 
the spacecraft spin period. 

MAGNETOMETERS 

Magnetometers, which define cone angles for use in attitude determination, measure the 
magnetic field components along three mutually orthogonal axes. The magnetic field com- 
ponents from such a triaxis magnetometer define the magnetic field vector in spacecraft 
coordinates for a particular universal time. The field vector in inertial coordinates is also 
known at that time from spacecraft orbit and magnetic field model computations. These 
two separate vector representations can be related to derive a unique cone angle between 
the inertial coordinate field vector and a possible spin axis orientation. 
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Figure B-1. Nine-bit Solar Sensor 




One triaxis magnetometer consists of three Shonstedt induction magnetometers, each 
aligned along a spacecraft axis with a sensitivity of ±27.853 milliamperes/meter (±350 milli- 
oersteds). The three field components must be sampled sequentially by the spacecraft 
telemetry, but the magnetic field vector changes slowly with respect to the telemetry 
sampling, such that no appreciable error is introduced. The time of measurement of the 
components is determined by the known time lag from sampling the magnetometers to 
the samples entering the telemetry, by the known spacecraft relative time of the telemetry 
magnetometer sample words, and by the cycling speed of the spacecraft telemetry related 
to universal time (deduced on the ground after transmission of the data). 

Spacecraft orbit positions for given universal times are computed from orbit elements pro- 
duced from ground observations of the spacecraft. These orbit positions have associated 
magnetic field vectors which are closely approximated by complex field models. The models 
are especially accurate for spacecraft altitudes between 300 and 2000 km. The field model 
vectors are transformed to inertial coordinates, and serve as the known reference vectors for 
cone-angle evaluation. 

The cone angles are constructed from the magnetic field components measured by the tri- 
axis magnetometer. If the components in spacecraft coordinates are represented by M x , 

M , and M , then the cone angle /x (figure B-2) is determined from 
y z 

M z 

cos n = — 

V^ x 2 + M y 2 + M z 2 

where, as before, the spacecraft z-axis is the spin axis. For a magnetic field vector Ft (in 
inertial coordinates), corresponding to the spacecraft external magnetic field at spacecraft 
position^., the spin axis must have coordinates of a point on a circle of radius /x about the 
vector it. (figure B-3). 



Figure B-2. Magnetic Field Vector Figure B-3. Magnetic Field Vector in Inertial Co- 

ordinates with Cone of Possible Axis Positions 
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EARTH HORIZON SENSORS 


One class of earth sensors used in cone-intersection attitude-determination systems consists 
of a telescope with a narrow field of view, mounted at a known angle y (~90°) to the 
spacecraft spin axis. It contains a sensing element to detect the change in visible light 
intensity between the sunlit earth and space or between the sunlit earth and the shadowed 
earth, and also the associated electronics to control the counters required for attitude solu- 
tions. Additional data, however, are required before the earth-horizon-sensor readouts lead 
to the cone angle of interest, the angle between the spin axis and the earth nadir. These 
data include (1) the earth center vector, (2) sun sensor data, (3) relations between the sun 
and the earth, (4) orbit information, and (5) reduction algorithms. The sun relations and 
the algorithms to find the desired relations are necessary because the earth is not a point 
source, as the sun can be considered to be, and the earth sensor readings are in relation to 
i points on the horizon, not to the earth center. 

The earth center vector is easily derived from orbit determination, since the position vector, 
of the spacecraft referenced from the earth center is opposite in sign to the earth center 
vector, li, referenced from the spacecraft. Thus, the components 

P x = -E x , P y = -E Y , P z = -E z . 

By this the earth center coordinates are known. Then, if the earth sensor, at a known 
angle y to the spin axis, produced a pulse only on sighting the center of the earth, the 
cone-intersection parameters would be known. Unfortunately, the earth sensor provides 
two pulses triggered by two unknown points on the earth horizon (figure B-4) or one 
unknown point on the earth horizon and one unknown point on the earth terminator, 
where the terminator is the division between the sunlit earth and the shadowed earth. 

Also, by means of the counters measures are made of the earth chord between the two 
unknown horizon/terminator points, of the sun-to-earth horizon time (earth time), of the 
telemetry-to-spacecraft sensor time (sun time), and of the successive sun sightings (the spin 
period). The earth chord is often inaccurately called “earth width” in the literature and 
that convention will be followed here. 

The telemetry readouts of the counter provide the means of relating the center-of-earth 
vector to the unknown horizon vectors. The counters are driven by a constant oscillator 
of C n counts/second. Thus the spin period, SP, is determined in seconds from the number 
of counts in the spin period counter, C sp , and the count rate, G n , by 

SP = C sp /C n . 

Similarly, the earth time is E t = C E t /C n , the earth width is EW = C Ew /C n , and the sun time 
is S ( = C St /C n , all in seconds of time. The earth time is expressed as the angle between the 
sun/spin axis plane and the earth horizon/spin axis plane through the relation 

E t 

v (in degrees) = — X 360° . 
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Figure B-4. Earth Horizon Sensor Relations 

Similarly, the earth width is written as an angle relation by 

Ew 

p (in degrees) = — 

SP X 360° 

The sun sensor produces the sun angle jS between the spin axis and the sun. Orbit informa- 
tion supplies the angle r\ between the sun and the earth center (negative spacecraft position 
vector). The angular radius of the earth p is derived from 

Re 

sinp = — 

R 

where R e = the radius of the earth, and 

R = the spacecraft position magnitude. 

NOTE: p is half the true angular earth width. 
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With these parameters, the earth nadir angle 5, between the spin axis and the earth center, 
can be determined as required for a second cone angle. Two earth-sun-spacecraft relations 
can exist, which allow and affect algorithms necessary to find 5. If from the spacecraft the 
earth is completely shadowed, no solution of 5 is possible; if the earth is completely sunlit, 
a comparatively simple algorithm making use of the half angle of n is used. If an earth 
terminator is visible, additional auxiliary angles between the earth horizon and the sun are 
formed, since the earth terminator is not directly useful. The two relations-for two hori- 
zons and one horizon/one terminator— are developed in two different algorithms in the 
section, “Some Non-least Squares Methods of Solution.” 

One other problem in dealing with earth-horizon sensors concerns the angular diameter of 
the telescope field of view. This angular diameter or scanning angle r causes the sensor 
photodiode pulse to trigger the earth horizon counters when a certain percent of the field 
of view is filled by the sunlit earth. In general, the pulse is not exactly coincident with the 
center of the field of view crossing the earth horizon. Therefore, the resultant counts are 
proportionately larger or smaller than an exact count. This generates a systematic error, 
which with some difficulty is largely eliminated. 
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APPENDIX C 

RESOLUTION OF SPIN-AXIS AMBIGUITY 


If a spacecraft has two sensors-a vector sensor and a Q vector sensor-to provide measure- 
ments for the determination of the spacecraft spin axis S, the vector relations of the section, 
“The Unique Solution, Three-constraint (PQV) Approach,” are applicable, that is, 

P • S = cos /3, Q • S = cos 5, P X Q = V = sin rj V, 

Q X P = -V, P • Q = cos 17 

where S has a two-solution ambiguity. In addition to these vector relations, 
period (SP) of the spacecraft, the time (Tj,^ Q ) between sensor sightings, and 
between the two sensors in the spin (xy) plane must be known or measured. 

S 2 be the two separate spin axis solution vectors, with 

V • Sj = -V • S 2 = cos a . 


the spin 
the angle (e) 
Let Sj and 


The convention is to define the spacecraft spin axis as the axis, viewed from above the 
spacecraft, about which the spacecraft spins ccw. The adjusted time is 


T p ^ q T' p ^ q + 36QO X SP 


where e/360° X SP is positive whenever the P vector sensor lags behind the Q vector sensor 
as the spacecraft spins. Typically, an attempt is made to align the sensors in the spacecraft 
xz plane to eliminate or minimize this factor. 

The selection of the correct spin axis solution is made by the following test: 


If a <90° and T p _> Q <HSP 
or 

if a > 90° and T p _ Q > 'A SP, 

then S = S x ; otherwise, S = S 2 . Refer to figure C-l for geometrical relations. 
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APPENDIX D 


THE METHOD OF LEAST SQUARES 

A set of n linear equations, 

measured calculated 

Y. « aXj + b (i = l,2,...,n), 

is overdetermined if n > 2, in that just two relations of measured Y. at a given X. are needed 
to solve unambiguously for a and b, but errors in the measured values Y. will not result in 
the same coefficients a and b for every choice of (X., Y.) (figure D-l). The requirement of 
a single set of the slope and Y intersection coefficients a and b, which makes the defined 
line closest in some sense to all of the measurements, is a problem often encountered. 

A least-squares best fit of the data is one widely used method to determine the coefficients, 
not only in the simple linear case, but also in many other applications. The linear case will 
motivate the discussion, but the results in matrix notation are general and can be applied 
wherever the matrix forms can be defined with a linear set of coefficients. 

The method called least squares refers to the mathematical procedure followed in which 
the differences or residuals (figure D-2) between the measured and calculated values, 

Y. - (aX. + b) = r ; (condition equations), 

are squared and summed over all points, the sum of the squares 


n 



i= 1 


is minimized with respect to the coefficients a and b, and the resulting two equations are 
solved for a and b. Since the sum of the squares is a function of a and b, 

n 

y; i \ = f o> b )> 

i= 1 


Y 



Figure D-1. Measured versus Calculated Values Figure D-2. Residuals between Measured and Cal- 

culated Values 
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minimizing the function is done by taking the partial derivatives with respect to a and with 
respect to b and setting them to zero : 

3f (a, b) _ Q df (a, b) q 


da 


9b 


The particular a and b for which these conditions hold simultaneously define the function 
f (a, b) such that f is a minimum (that is, f has a greater value for any other choice of a and 
b). 


Carrying out the partial differentiations and solving for a and b, 
9f (a, b) 

= L 

i= 1 i= 1 


9a 


n n 

= 2 ^ (Y. - b - aXj) (-X.) = 0 = 2 ^ (aX 2 + bX. - X.Y.) 


or 


9f (a, b) 
9b 


n ii 

= 2 ^ (Y i - b - aX ; )(-l) = 0 = 2^ (aX. + b - Y.) 
i= 1 i= 1 


i=l i=l i= 1 


n 


n 


■E x ' tnb ’ Y y > 

i= 1 i= 1 


(normal equations) 


wherefore 


n n n 

n L Y ' x ‘- E y ‘ L x ' 


i= 1 


i= 1 


i= 1 


a = 


n it 

X> 

i= 1 i- 1 


n n n n 

Y x ? Y y , - Y x . y . Y x < 


b = 


i= 1 


i= 1 


i= 1 


"Y x '- Y x ‘ 


i= 1 


i= 1 
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The same results are produced by an equivalent but more convenient matrix least-squares 
approach. The initial set of linear equations is given by 



where the approximations and residuals do not enter explicitly into the matrix approach.* 
In representational form, 

(J !)'(»? 2) (2 XI) (C ° ndl,km ' eq “ i<, “)' 

Premultiply by the transpose of X to square the prematrix of C 
X T Y = X T X C 

. gXnXnXl) (2 X n) (n X 2) (2 X 1) („ orma | equatiom) . 

(2X1) | (2 X 1) 

If the (2 X 2) X T X matrix is nonsingular, it has an inverse and 

(X T X) 1 X t Y = C 
(2X2) (2 X 1) (2 X 1) 


which in expanded form is 



♦Smith, G. A., The Theory and Application of Least Squares, GSFC TMX 63127, December 1967. 


D-3 



This is the result obtained by squaring residuals, summing, minimizing, and solving for a 
and b. 


Certain equations involving nonlinear parameters can be solved for these parameters by the 
least-squares method using an iterative process. Beginning with a set of nonlinear equations 
and an initial estimate of the desired parameters, the equations are linearized and solved for 
corrections to the initial estimate. The corrections are added to make a new estimate and 
the procedure is repeated if the corrections are not less than preselected limits. 


Since the required parameters— coefficients such as a and b or components of other rela- 
tions-are not involved in a simple, linear fashion in the defining equations, these defining 
equations are linearized by expanding in a Taylor’s series about the required parameters, 
truncating the series after the linear terms, and evaluating the function at the estimated 
solution. Assume a general function, f. = f. (Y., Xj, A k ) = 0 (i = 1, 2, . . ., n; k = 1 , 2, . . . , 
m), is composed of a set of n equations which are nonlinear in the m coefficients A k . 
Expanding in a Taylor’s series about A k evaluated at the estimate A° , 


f,(Y 1 .X 1 .A k ) = f 1 (Y 1 ,X 1 ,A” t )* 


L 

lr = 1 



I AA k + 0(AA k ) 2 
A = A° 

A k 


+ . . . 


where 


AAfc " (A„ - A° k ) 


For additional notational simplicity, let 

f.(Y.,X.,A° k ) = f" and 


/dl \ df. 

vv 0 = "k 


A = A 
A k A k 


Then the expansion truncated after the AA k terms is 


m 

f,<Y 1 .X 1 .A k )*fr ♦ 


k= 1 


Switching to matrix notation to complete the solution of the corrections AA k , 
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(nX 1) (nXm) (mXl) 

or 

-F = PA 

where F = the matrix of function equations evaluated with the estimated coefficient 
A° 

A k> 

P = the matrix of partial differentiations, and 

A = the matrix of coefficient correction elements. 

The procedure for solving for A is identical to the first case of solving for C since the equa- 
tions are now linear in A, such that 

A = -(P 1 ?) 1 P t F . 

The correction elements are compared with the preselected limits E k and depending on the 
comparison results, different actions are taken. If some AA k > E k , then new estimates 

A °k = A °k + AA k 

are formed and the process is repeated. If all A k < E k , then the final, acceptable solutions 
are formed 

A k = A k + AA k • 


NASA-Langley, 1975 
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